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ABSTRACT 


This thesis explores Benders decomposition for solving interdiction 
problems on electric power grids, with applications to analyzing the vulnerability 
of such grids to terrorist attacks. We refine and extend some existing 
optimization models and algorithms and demonstrate the value of our techniques 
using standard reliability test networks from IEEE. 

Our implementation of Benders decomposition optimally solves any 
problem instance, in theory. However, run times increase as Benders’ cuts are 
added to the master problem, and this has prompted additional research to 
increase the decomposition’s efficiency. We demonstrate empirical speed ups 
by dropping slack cuts, solving a relaxed master problem in some iterations, and 
using integer but not necessarily optimal master-problem solutions. These mixed 
strategies drastically reduce computation times. For example, in one test case, 
we reduce the optimality gap, and the time that it takes to achieve this gap, from 
16% in 75 hours to 5% in 16 minutes. 
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DISCLAIMER 


The reader is cautioned that computer programs developed in this 
research may not have been exercised for all cases of interest. While every 
effort has been made to ensure that the programs are free of computational and 
logic errors, they cannot be considered fully validated. Any application of these 
programs without the additional verification is at the risk of the planner. 
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EXECUTIVE SUMMARY 


This thesis extends existing optimization models and methods for 
analyzing the vulnerability of electric power grids to terrorist attacks. 

Electric power systems are critical to the United States’ economy and 
security. After the events of September 11, 2001, and the United States’ ensuing 
war against terrorism, the fear of retaliation against critical infrastructures has 
become a major concern for security analysts. The vulnerability of the electric 
power systems to physical disruptions heads the list of concerns, because the 
U.S. transmission grid has not expanded as quickly as demand has over the last 
decade, and thus the system has become “brittle.” This brittleness (vulnerability) 
became more evident on 14 August 2003 when a relatively modest number of 
system malfunctions caused a blackout of the northeastern region of the country, 
and raises this question: How much worse might the blackout have been had it 
resulted from a well-planned terrorist attack? 

This thesis first considers certain modeling issues of the problem of 
optimal interdiction, and then focuses on techniques for solving the proposed 
models. This should help us identify the most critical components of the U.S. 
power grid, and ultimately help determine effective protective measures to 
improve the system’s robustness. 

We refine and extend the existing optimization models and algorithms of 
Salmeron et al. that identify critical system components (e.g., transmission lines) 
by creating maximally disruptive attack plans for terrorists, who are assumed to 
have limited offensive resources. Most importantly, this thesis exploits bi-level 
structures embedded in the interdiction models to enable faster solutions through 
the use of Benders decomposition. 

We first validate the DC power-flow model that underlies in the interdiction 
models. By manipulating a few equations, the inherent nonlinearities in these 
models are substituted by linear forms, converting those models into (linear) 
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mixed-integer problems (MIPs). We explore Benders decomposition for solving 
the MIPs, and successfully apply it to two standard reliability test networks from 
IEEE. 

Although we solve the IEEE test cases to optimality, we notice that the 
efficiency of the Benders’ master problem deteriorates as (a) the number of 
components in the test cases increases, and (b) the number of Benders’ cuts in 
the master problem grows. We investigate several techniques in order to 
improve the algorithm’s speed. First, we demonstrate faster convergence by 
solving the mixed-integer master problem exactly every /c th iteration only, and 
solving its easier, linear-programming relaxation otherwise. Importantly, the 
algorithm thus modified always maintains a valid upper (optimistic) bound, just as 
the original does. Since relaxed (non-integer) master-problem solutions cannot 
be used by the subproblems to compute lower (pessimistic) bounds, we exploit 
sub-optimal integer solutions to the master problem every / 77 th iteration to improve 
the lower bound. In one of our two test networks, this combined techniques 
reduce computational time by one third. We also show that limiting the number 
of cuts in the master problem can speed convergence, and study several 
strategies for keeping or deleting cuts. Results are mixed, but improvements can 
be dramatic. For example, in the second of our test networks, by limiting the 
number of cuts to 500 and solving a relaxed master problem nine of every ten 
iterations, we achieve an optimality gap under 5% in 16 minutes; this compares 
to a 16% gap in 75 hours using the original decomposition algorithm 



GLOSSARY 


The following briefly defines the general concepts related to electric power 
networks and interdiction. Some of these definitions are from (or strongly based 
on) the glossaries provided by Energy Information Administration (EIA 2003), 
Elec-Saver (2003) and SIEMENS(2004). 

AC Power: Power associated with alternating current circuits. 

Admittance: Property that allows the flow of electrical current through 
reactive circuit elements under the action of a potential difference. Admittance is 
the reciprocal of impedance. Admittance equations establish the relationship 
between current, impedance and voltage. 

Alternating Current: Current that periodically reverses direction. 

Bus (or Busbar): A heavy, rigid electrical conductor that makes a common 
connection between several electrical circuits. 

Capacitance: The property of a circuit that allows it to store an electrical 
charge. 

Case: A set of data to be analyzed, along with the results of the analysis. 
Data consists of power grid components (lines, buses, generators, substations, 
etc), physical data (i.e., such as line impedances, generating capacities, etc.), 
non-physical data (e.g., interdiction resource, optimization parameters, etc.). 
Results include the interdicted components, associated generation and power 
flows, and load shedding for every period, among others. 

Costumer Sector: A type of load with specific requirements (e.g., amount 
of power demand and cost for failing to provide it). 

Current: The flow of electrons in a circuit. Current is measured in 
amperes. 

DC Power: Power associated with direct current circuits. 
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Direct Current: Current that flows in one direction 


Disruption: The cost of power shed (in dollars per hour) or the cost of 
energy shed (in dollars), as a consequence of interdiction. 

Energy: The result of integrating power over time. 

Energy Shed: Amount of energy that cannot be supplied to the load (one 
or several customer sectors ) over the course of a given period. 

Generation (of electricity): The production of electric energy through the 
transformation of other forms of energy. The amount of electric energy 
produced. 

Generating Unit (or Generator): Any combination of physically connected 
generator(s), reactor(s), boiler(s), combustion turbine(s), or other prime mover(s) 
operated together to produce electric power. 

Impedance: The total opposition to alternating current. Impedance is the 
vector sum of resistance and reactance. The unit for impedance is the ohm. 

Inductance: The property of an electrical circuit that causes it to oppose 
changes in current. 

Interdiction Plan: Specific subset of the electric system equipment that 
might be interdicted by terrorists. Optimal or near-optimal interdiction plans are 
identified in for given interdiction-resource scenarios. 

Interdiction Resource: A numerical value associated with a mathematical 
expression that represents the capacity of terrorists to carry out attacks. For 
example if such an expression has the form: “(3 x total number of attacks to 
buses) + (1 x total number of attacks to lines) < 5,” then “5” is the interdiction 
resource. 

Line: See “Transmission Line.” 

Load: Demand for electric power, measured in watts, at a specific point in 

time. 
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One-line Diagram: Schematic drawing of an electrical power system that 
uses graphical symbols to represent electrical equipment such as buses, 
generators, loads, transmission lines and transformers. It may incorporate 
numerical values for the system, such as line power flows, generating unit 
outputs, bus voltages, etc. 

Phase angle: The angle by which the sine curve of the voltage in a circuit 
element (or a combination of elements) leads or lags the sine curve of the current 
in that circuit element(s). 

Period (of restoration): Each of the stages that an electric power network 
undergoes following an attack, as interdicted components are repaired or 
replaced over time. 

Power (transmission): The transfer rate for energy. Unit of measurement is 
watts. Power is sometimes used loosely to refer to electricity. 

Power Shed: Amount of power that cannot be supplied to a load or loads 
(one or several customer sectors) at a specific point in time. 

Reactance: Opposition to the flow of alternating current. Measured in 

ohms. 

Reactive Power: Power associated with inductance or capacitance. 

Resistance: The propensity of a circuit to oppose current flow. It is 
measured in ohms. 

Scenario: A particular value of the interdiction resource. In our 
formulation, we analyze the worst-possible Interdiction Plans for a given 
scenario. 

Slack Bus: (Also called swing bus.) A bus whose (initially unspecified) 
power generation must be determined in order to match supply and demand in 
the network. 

Substation: Facility with equipment that switches, changes, or regulates 
electric voltage and current. 


XXI 



Transformer: A static electrical device which, by electromagnetic 

induction, regenerates AC power from one circuit into another and/or changes 
the voltage of alternating current. 

Transmission: The movement or transfer of electric energy over an 
interconnected group of lines and associated equipment between points of 
supply and points at which it is transformed for delivery to consumers, or is 
delivered to other electric systems. Transmission is considered to end when the 
energy is transformed for distribution to the consumer. 

Transmission Line (High-voltage): The high-voltage (> 69 kV in the U.S.) 
conductors used to carry electrical energy from one location to another. 

Transmission System (also referred to as Electric Power Grid and Electric 
Network in this thesis): An interconnected group of electric transmission lines 
and associated equipment for moving or transferring electric energy in bulk 
between points of supply and points at which it is transformed for delivery over a 
lower-voltage distribution system to consumers, or is delivered to other electric 
systems. 



I. INTRODUCTION 


This research continues the development of mathematical models and 
optimization methods for planning electrical power grids that are robust to 
terrorist attacks. We focus on solving an interdiction model, which is represented 
as a mixed-integer program (MIP), through the use of Benders decomposition. 


A. THE ELECTRIC POWER GRID IN THE U.S. 

The electric power grid in the United States (U.S.) consists of over 10,000 
generating units with a total production in excess of 760 giga-watts (GW), and 
over 700,000 miles of transmission lines, all controlled by approximately 100 
control centers. (Actually, only modest interconnection capacity exists between 
three main sub-systems in this grid, namely the Eastern interconnection, the 
Western interconnection and the Texas interconnection. Therefore, “the U.S. 
power grid” may be viewed as three, nearly independent systems.) 

In the past, utilities have been concerned with the grid’s vulnerabilities to 
isolated natural disasters and unscheduled, technical outages caused by 
equipment failures. To date, protecting the electric grid against multiple, nearly 
simultaneous failures has not been a high priority for utilities because of the 
expense involved and the relatively low frequency of such events. However, in 
the current era, a set of nearly simultaneous failures might be the precise 
objective of a terrorist attack. “Low frequency” may no longer be so slow. 

As a result of the threat that terrorism currently poses to critical 
infrastructure of all types in the U.S., understanding the physical vulnerability of 
the electric grid has become more important not only for the electric power 
industry, but also for national security (NSTAC 1997). In fact, a National Security 
Assessment (1997) states that man-made physical destruction is the greatest 
threat facing the electric power infrastructure. 

The U.S. electric power grid is plagued by a number of inherent 
weaknesses. Among the most significant of these weaknesses are the traversals 
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of sparsely populated areas, which readily expose the system to malicious 
attacks, and the decreasing reserve in transmission capacity. 

The U.S. transmission network is built with some redundancy and multiple 
protective devices; however, current reserve levels in transmission capacity may 
be insufficient to back up multiple failures in some key components, as 
exemplified by the 14 August 2003 power-system failure (which cost our national 
economy between $7 and $10 billion dollars; ICF Consulting 2003). Not 
surprisingly, this failure brings additional speculation about terrorist threats to the 
U.S. power grid (Germain 2003), and raises important questions about the 
resilience of the U.S. power grid and the possibility of widespread loss of power 
to customers. 

The August 14 th blackout could be financially insignificant when compared 
to the effects of a coordinated terrorist attack on several key facilities during a 
period of peak load. ICF Consulting (2003), which analyzed the cost of the 
August 14 th blackout, states that a terror-induced blackout could prove 
significantly more costly and have debilitating impacts on the affected region as 
well as the entire country. Additionally, a simulated terrorist attack on the Pacific 
Northwest's power grid was recently conducted under the project name “Blue 
Cascades.” An analysis of that simulated attack showed that a real attack, if 
successful, could wreak havoc on the nation's economy, shutting down power 
and productivity in a domino effect that would last for weeks (Allen 2003). 

The discussion above motivates the development of optimization models 
to represent the problem of terrorists attacking a power grid. By studying how to 
attack power grids, one can gain insight on how to protect them. 

This research follows Salmeron et al. (2003, 2004), who develop bi-level 
mathematical models to identify maximally disruptive attacks for terrorists, who 
are assumed to have full access to power-system data, but limited offensive 
resources. Salmeron et al. initially introduce a static interdiction model and then 
extend it to capture the dynamics of system operation as the grid is repaired after 
an attack (which can make slow-to-repair grid components more attractive 
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targets). An attack is optimal in the first model if it maximizes the total, 
instantaneous, unserved demand for power; in the extended model, it is optimal if 
it maximizes the total cost of unserved demand for energy. These models have 
been tested using standard reliability test systems (RTSs) drawn from IEEE 
(1999-1, 1999-11). 

This thesis makes no attempt to extend the modeling approach for power 
flows used by Salmeron et al. Thus, some values representing power-system 
behavior (such as reactive power flows, line losses, voltage magnitudes, 
transformer tap positions, etc.) are still disregarded and/or are assumed fixed in 
this work. However, we assess the accuracy of that approach by comparing 
results from our approximating “DC power-flow model” with those provided by a 
full AC power-flow model. We use the AC model embedded in the PowerWorld 
Simulator (PowerWorld 2003) to carry out this comparison. 

Salmeron et al. (2003) propose heuristic algorithms for finding an 
approximate solution to their interdiction models. These heuristic solutions are 
generally within 10% of optimality but in some instances that gap can reach 25%. 
The possibility of large gaps was discovered after publication of the referenced 
paper, when the authors were able to solve equivalent Mixed Integer Linear 
Program (MIP) formulations exactly (Salmeron et al. 2004). These MIP 
formulations have been reviewed and consolidated in the framework of this 
thesis 

In theory, a MIP can be solved exactly by standard techniques such as 
“branch-and-bound” (e.g., Wolsey 1998, pp. 95-107). However, these 
techniques cannot always solve the large-scale problems found in practice, and 
this motivates the search for alternative techniques with better scalability. 
Benders decomposition (Benders 1962) is a well-known method for solving large- 
scale MIP models, and this thesis investigates its use to solve the interdiction 
problem efficiently. Previous application of Benders decomposition to interdiction 
problems can be found in Cormican (1995) and Israeli and Wood (2002). 


3 



B. THESIS OUTLINE 

Subsequent chapters in this thesis are organized as follows: Chapter II 
introduces the DC power-flow model that is the basis for the interdiction model. It 
also defines the interdiction model and presents the heuristic and MIP 
approaches developed previous to this work to solve the problem. Finally, 
Chapter II introduces Benders decomposition and its associated algorithm in the 
framework of our interdiction problem. Chapter III validates the DC power-flow 
model through comparison with an AC power-flow model. We compare power 
flows for different interdiction plans on IEEE’s “One Area RTS.” Chapter IV 
details the application of Benders decomposition to our interdiction model. A 
three-bus example demonstrates the decomposition process, which is then 
applied to the interdiction problem without system restoration. Computational 
results are presented. Chapter V presents the Benders decomposition for the 
interdiction model with system restoration over time. Chapter VI discusses 
refinements to the algorithm and provides computational results. Chapter VII 
summarizes results and recommends topics for future work. Appendices A and 
B show the full derivation of the interdiction models, which are introduced in 
Chapter II and further described in Chapters IV and V. Appendix C describes the 
linearization of cross products for the three-bus problem developed in Chapter 
IV. 
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II. PRELIMINARIES 


This chapter introduces the DC Optimal Power Flow model (DCOPF) that 
is the basis for the development of two interdiction models: without system 
restoration (see Section B), and including system restoration (see Section C). 
Both Sections B and C describe heuristic and MIP approaches explored prior to 
this research for solving the interdiction models. Section D introduces Benders 
decomposition in the context of our models. 

A. INTRODUCTION TO DCOPF 

At the heart of the interdiction models is a DC power-flow model (DCPF), 
so we first provide some background on this and related models. DCPF 
simplifies the so-called “full AC power-flow model” (ACPF), which is generally 
accepted as a valid representation of power flows under steady-state conditions 
and symmetry (Frauendorfer et al. 1993). The DCPF representation entails 
various assumptions which may be acceptable in the context of security analysis 
(Wood and Wollenberg 1996). For example, DCPF disregards reactive power 
effects and assumes nominal voltage magnitudes at the buses. 

DCPF can be used to construct an optimization problem, DCOPF, which 
optimizes a merit function subject to DCPF constraints. The merit function 
typically measures generating costs. In contrast, our DCOPF not only minimizes 
generation costs, but also the penalty associated with unmet demand (“load 
shed”), because we cannot guarantee that all demand will be met in the 
presence of interdiction. The DCOPF formulation from Salmeron et al. (2003) is 
presented below for completeness. (Definitions of terms can be found in the 
glossary.) 

Sets: 

iel, set of buses 

ge G , set of generating units 
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le L, set of transmission lines 

ceC, set of consumer sectors 

se S, set of substations 

iel s , subset of buses at substation 5 

ge G l , subset of generating units connected to bus i 

le Lf s , subset of lines connected to bus i 

leL s s ub , subset of lines connected to substation 5 (including 
transformers, which are represented by lines) 

G*cG, I’ci, I* e / , S* c5 , interdictable generators, lines, buses, and 
substations, respectively. These are “interdictable 
components.” 

Parameters (units, if applicable): 

oil), d(l), origin and destination buses of line / respectively (more than 
one line with the same o(l ), d(l ) may exist). 

Kg. ), bus for generatorg, i.e., ge G i(g) 

s(i), substation se S associated with bus ie I s 

d ic , load of consumer sector c at bus i (MW) 

Pj Line , transmission capacity for line / (MW) 

P g en , maximum output from generator g (MW) 

r n x, , resistance, and reactance of line / respectively (Q.). (We 

assume x, » r, ) 

B,, series susceptance, calculated as B, =x,/(rf +xf) (j/^) 
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generation cost for unit g ($/MWh) 


f ic , load-shedding cost for customer sector c at bus i ($/MWh) 


Decision variables (units): 

P ° en , generation from unit g (MW) 

P , Line , power flow on line / (MW) 

S ic , load shed (unmet) for customer sector c at bus i (MW) 

0 j , phase angle at bus i (radians) 


Remark: All units are converted to per unit (p.u.) values for a base load of 
100 MW. (Remark: We use this same value in all implemented models.) For 
example, to convert r, to per-unit values, we consider the transmission line (/) 

nominal voltage rating E, (in kV). Then: 


r, (per unit) = /} (ohms) * 100 (MW)/(line rating(kV)) 2 = 


100/; 

( MVA a 

100/} 

f mXa _q ' 

100/} | 

f a 

Ef 

l (kV) 2 J 

" Ef 

l (X 2 yV J 

' Ef 1 

l v J 


100r 

—T^-p.u. (since 1V=1 Ax1£2) 


The same process is carried out for reactance (x ,) conversions. For the 

rest of MW magnitudes, we use per-unit values after dividing by the 100 MW 
base load. 

The formulation of DCOPF is: 


(DCOPF): min 

pGen pLine ^ q 


SVf+SS/A 

g i c 


s.t. 


P l Line =B l (0 oa) -0 d{l) ) 

V/€ L 


(2.1) 

^ ' pGen ''y ' pLine _|_ 

z pr=Z(d„-s iI ) 

Vz 

(2.2) 

g l\o(l)=i 

l\d(l)=i c 



pLine ^ pLine ^ pLine 

V/eI 


(2.3) 
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V/,Vge G t 
Vz',c 


(2.4) 

(2.5) 


pGen ^ pGen 
—8 ~ g 

0 <S,<d,„ 


<P 


Gen 


This model assumes fixed (nominal) voltage magnitudes and does not 
represent shunts, reactive power flows, power losses, DC lines, transformer tap 
positions, phase transformers, and other features that could destroy the linear 
structure of the approximation. Discrepancies in power flows in DCOPF and a 
full ACPF are discussed in Chapter III. 


B. THE INTERDICTION MODEL WITHOUT SYSTEM RESTORATION 


1. Interdiction Model 

After some manipulation, DCOPF above can be stated in standard form 
as: 


(DCOPF): min cy 

y 

s.t. Ay = b (2.6) 

v > 0 

where y represents generation outputs, load shedding, power flows and phase 
angles variables, along with slack variables associated with constraints (2.3)- 
(2.5). 


The interdiction model extends (2.6) in order to choose the most 
disruptive, resource-constrained interdiction plan, Se A, where A is a set of 
binary vectors representing possible terrorist attack plans. An interdiction plan is 
represented by the binary vector 8 , where 8 e is 1 if component e of the power 

grid is attacked and is 0 otherwise. We assume that terrorists attempt to 
maximize the minimum post-interdiction “disruption” which is measured as the 
generating cost plus the load-shedding cost (penalty) evaluated by DCOPF. 


A simple representation of the max-min Interdiction of the DCOPF (I- 
DCOPF) model is shown below: 
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(I-DCOPF) : max min cy 

Se A y 

s*t. g(8,y) = 0 (2.7) 

v>0 

where g(S,y) represents a set of constraints, some of which involve functions 
that are nonlinear in ( 8,y ) (see Appendix A.1). For a given interdiction plan, 8 , 
the inner minimization is still a DCOPF model which is, of course, linear in y. 

Prior to this thesis, model (2.7) was tackled by (a) using a heuristic 
algorithm and by (b) reformulating (2.7) as a MIP which was solved with a 
standard solver. This thesis contributes to (b) by refining the MIP formulation 
that serves as basis for the Benders decomposition approach. In order to place 
this thesis’ contributions in context, the following two subsections summarize 
approaches (a) and (b). 

2. Heuristic Approach 

A heuristic algorithm introduced by Salmeron et al. (2003) solves model 
(2.7) approximately, and has been tested using small- and medium-scale 
networks drawn from the IEEE Reliability Test Data (1999-1, 1999-11). This 
heuristic solves DCOPF for a given grid configuration (i.e., for a given interdiction 
plan). Then it assigns “values” to interdictable electrical components in the grid 
(lines, buses, substations, etc.) based on present and previous flow patterns 
associated with the components. Finally, the heuristic maximizes the value of 
the components to be interdicted while excluding previously explored solutions. 
This process is repeated for a pre-specified number of iterations. This 
decomposition-based approach has been empirically proven to obtain good, but 
not necessarily optimal, interdiction plans. Figure 1 depicts the basic cycle of the 
algorithm. Additional details on the heuristic can be found in the cited reference. 
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Figure 1. Framework for a heuristic interdiction algorithm. (From Salmeron et al. 

2003.) 

The inability to prove optimality of the heuristic solution (unless all possible 
solutions are enumerated explicitly) has led to the search for alternative 
representations of (2.7) that are amenable to solution by MIP techniques, such 
as branch-and-bound, as described in the following section. 

3. The MIP Approach 

a. Overview 

The development of a MIP model enables validation of the heuristic 
solution and the development of exact decomposition algorithms that are 
frequently used solve large-scale problems. (Typically the solutions are 
approximate, but then large-scale systems are almost never solved exactly by 
any technique.) The derivation of the MIP begins with the linearization of 
constraints in g(S,y) = 0, followed by the dualization of the inner minimization 
problem in (2.7). This procedure converts the max-min problem into a nonlinear 
maximization problem which is linearized through additional steps. Both of these 
linearization techniques are fully explained in Salmeron et al. (2004), but are 
briefly described in the following for continuity. 

b. Derivation 

The nonlinearities in the expression of g(S,y) = 0 in model (2.7) are 
associated with admittance equations in the presence of interdiction variables £. 
These have the form: 
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P, = (1 - 4" XI - <0(1 - 8m > 


no-o 

V*i ,e 4 y 


B,(0 oU) -e dO \VleL (2.8) 


Note that in the interdiction model, the original admittance equation 
(DC.1) for every interdictable line / is modified as (2.8) to account for potential 
interdiction (directly or indirectly) of the line, in which case the admittance 
equation should not be enforced. This leads to nonlinearities in the form of 
multiple cross-products. 


Following Salmeron et al. (2004), (2.8) can be fully linearized. For 
example, one can convert the nonlinear constraint P = B(0 a -0 b )(l-S l )(l-8 2 ) into 
the following two linear inequalities: 

P~B{O a -O b )< M(S l + 8 2 ) 

P-B(O a -9 b ) > -M(8 l +S 2 ) 

Here, M = P + Bd , where 6 is an upper bound on the maximum 
phase angle difference between buses a and b. (2.9) enforces P = B{6 a -G b ) 

when all the related £-variables are 0, and eliminates the constraint when any of 
these S variables is 1. The resulting model is still called l-DCOPF in this thesis, 
and is specified in detail in Appendix A.1. 


We proceed by taking the dual of the inner minimization in (2.7), 
assuming constraints in g(8,y) = 0 have been linearized as described above. 
This converts the max-min problem into a simple maximization. The resulting 
model is called DI-DCOPF (“Dual of l-DCOPF,” although only the inner 
minimization is dualized). The full formulation is given in Appendix A.2. This 
model’s objective function consists of a linear function of continuous vector n, 
along with nonlinear costs involving cross-products of type 8k. The new 
variables k correspond to duals of each of the balance and bound constraints in 
l-DCOPF. The terms associated with the 8k products arise from the line and 
generator capacity constraints. 
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The resulting problem can be stated, in brief, as: 

(DI-DCOPF): max max b*7T + b'v 

Se A n 

s.t. A n<c (v) (2.10) 

VG V(8,7r) 

Here, b* and b" are appropriately dimensioned coefficient vectors, 
7i is a vector of dual variables to constraints in l-DCOPF and v is a vector where 
each component is the result of a cross product of a variable 8 and a variable 
n. Variables in parentheses, here and elsewhere, are dual variables for the 
associated constraints; in this case the variables are y . 

To linearize a cross-product involving a 8- variable and a 
^-variable, we create an intermediate variable, e.g., v jk = 8 j n k , with the 

property that, 8 f = ()=> v /k = 0, and 8 t = 1 =» v jk = n k . Since the sign of n k is 
known (assume n k > 0 for demonstration purposes), we add the following linear 
constraints, represented as vgV(8,tt), to ensure the above relationships 
between 8. , 7t k , and v jk hold: 

v jk ^ x k Sj 

v ]k ^ K k 

v ]k = 8fa ~ ' v jk - n k ~(1 - 8j) ( 2 . 11 ) 

0<7T k < n k 

Vjk^ 0, 

where n k is an appropriate upper bound on n k . 

Let us drop the (j,k ) subindices for simplicity. Then, the following 
scheme depicts the validity of the representation in (2.11): 


12 



V<7rS 

V<71 

V>K-7t(\-S) 
0<7T<7T 
V > 0 


Remark: When n < 0, the linearization is carried out as: 

v jk ^ ~n k 8j 
v jk ^ x k 

v jk = 8pc k = < V jk < 7T k + K k (1 - 8)) (2.13) 

-71 k <7l k <{) 

v jk Z 0 

Finally, generation phase angles, power flows and load shedding 
are jointly denoted by the dual variable y . 

The dualization and linearization described above allow us to 
represent the interdiction problem as a MIP. The resulting model is called LDI- 
DCOPF (Linearization of DI-DCOPF; see Appendix A.3). 

C. INTERDICTION MODEL WITH SYSTEM RESTORATION 

1. System Restoration 

The model described in Section B is a static representation of power 
disruptions at a given point in time. That model does not represent the 
consequences of the variability in repair times of damaged system components 
and the increased cost of load shedding resulting from repair-time delays. 
Unless the outage duration of each interdictable component is the same, the 
interdiction model presented in Section B cannot capture this effect. For 




13 



example, Salmeron et al. (2003) show that if the interdictor’s decisions are driven 
by the goal of maximizing short-term power shedding (or its hourly cost), the 
resulting interdiction plan might be far from optimal had his or her goal been 
maximizing total energy shed (or its cost) until the system repair is completed. 
This issue is especially important when transformers and other difficult-to-replace 
equipment are subject to interdiction. 

Figure 2 depicts of the ideas above using two graphs. 



Figure 2. Power disruption and energy disruption. The model without restoration 
provides the optimal interdiction plan according to instantaneous picture of 
power shed (left). The model with restoration over time accounts for 
energy disruption (right). (From Salmeron et al. 2004.) 


2. Interdiction Model 

Salmeron et al. (2004) show how l-DCOPF (2.7) can be extended to 
handle the time associated with repairs. This resulting model seeks to maximize 
total cost of energy shed, accounting for the fact that the system will be 
progressively restored over time. The interdiction model with system restoration 
can be stated in brief as: 
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max min y,D l c v, 

& A V ^ 

s.t. gl (8,y t ) = 0, \/te T (2.14) 

y t >0, VteT 

This model uses interdiction constructs to couple instances of DCOPF, 
one for each system state that represents a stage or “time-period” of system 
repair: Note the time subindex t , where T is the set of time periods based on 
repair times of all components in the system (see Appendix B.1). That index is 
added to all of the inner decision variables to account for power flows in the grid 
and other decision variables that change over time. D t represents the duration 

of time period t, y t is the same y described in (2.6) but for every time period t, 
and g t is the resulting set of nonlinear equations in (S,y t ). Note the time- 
dependent constraints g,(8,y,) = 0 not only ensure that some components are 
out of service right after interdiction, but also guarantee that components will be 
in service again after their posted repair time. Thus, the constraints g,(8,y l ) = 0 
for different periods t are coupled by the variables 8. The full model derivation, 
described in detail in Salmeron et al. (2004), consists of the inner DCOPF model 
with system restoration (referred to as DCOPF-R) and the associated interdiction 
model (referred to as l-DCOPF-R); see Appendix B.1 for additional definitions 
and notation, and Appendix B.2 for the formulation. 

As in the case without system restoration, a solution to the interdiction 
model with system restoration can be obtained through a heuristic algorithm or 
through a MIP reformulation of the problem, as described in the following 
paragraphs. 

3. Heuristic Approach 
a. Description 

The Heuristic for a model with system restoration over time 
(schematically depicted in Figure 3) follows the approach of that for a model 
without system restoration. 
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For a fixed interdiction vector 8 (i.e., 8 = 8), the inner DCOPF-R 
subproblem provides the joint flow patterns for a number of system stages 
(restoration periods). The subproblem decomposes into |7j sub-subproblems 
(because when 8 is fixed, the subproblem is no longer coupled by time periods). 
Each subproblem consists of an instance of DCOPF with some subset of system 
components being “out of service.” Components that are out of service are 
determined “on the fly” as the algorithm solves a DCOPF model for every time 
period. 

Salmeron et al. (2004) redefine the concept of a component’s value 
in this dynamic model by factoring the time it would be out of service if 
interdicted. 


r 


Solve the DC-OPF-R for the given 
Interdiction Plan: 

Solve \T\ DC-OPF problems 


Based on present and previous flow 
patterns, assign an (energy-based) 
“Value” to each interdictable asset 


"S 


MP-R: Maximize the Value of the 
Assets to be Interdicted (excluding 
previously explored solutions) 


Figure 3. Interdiction algorithm framework with restoration. DCOPF-R decomposes 
into |7j sub-subproblems, each being an instance of DCOPF with some 
components “out of service.” (From Salmeron et al. (2004).) 


4. The MIP Model 
a. Overview 

The structure of l-DCOPF-R allows us to convert it into a standard 
MIP through a process of dualization and linearization, just as in the case without 
system restoration. 
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b. Derivation 

As the new interdiction model is a variant of the model without 
system restoration, where the inner model has a similar structure to a power flow 
model (replicated over time), the development of the MIP formulation is 
performed as in Section B: First, admittance equations are linearized (Appendix 
B.2). Next, the inner minimization problem is dualized (Appendix B.3). Finally, 
cross-products in the objective function are linearized. The resulting MIP, called 
LDI-DCOPF-R, is specified in Appendix B.4. 

D. BENDERS DECOMPOSITION 

Benders Decomposition is a well-known technique for solving MIPs 
(Benders 1962). This technique has proven effective in solving large-scale 
optimization problems (e.g., facility location problems, two-stage stochastic 
optimization problems, etc.) including interdiction problems (Cormican 1995, 
Israeli and Wood 2002). Benders decomposition is well suited when the 
assignment of a subset of variables, so-called “complicating variables,” yields 
easily solvable, disconnected and convex subproblems, such as the network-flow 
problems in Geoffrion and Graves (1974) and Brown and McBride (1984); these 
MIPs could probably not have been solved directly, i.e., using branch and bound. 
Based on the successful application of Benders decomposition to these similar 
problems, we implement it to solve LDI-DCOPF and LDI-DCOPF-R. 

Since the development of Benders decomposition is similar for both 
models, we focus on LDI-DCOPF. The MIP formulation of this model embeds a 
bi-level structure that makes Benders decomposition suitable for solving the 
problem: The first level is the interdictors’ level, which is modeled through binary 
decision variables. The second level is an instance of DCOPF (actually, its dual), 
where some components of the original power grid have been removed through 
the first-level decision variables. Because DCOPF is linear, it satisfies Benders’ 
theoretical requirement for convexity in the subproblem. 
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Our interdiction problem (LDI-DCOPF) may be stated, in brief, as: 

(P): z* = max c'v + c n 7l 

8 ,V ,7T 

s.t. A v v + A*x + BS = b (v) 

7T> 0 (2.15) 

v > 0 
Se A 

Note: This problem is referred to as LDI-DCOPF elsewhere in this thesis, 
but for simplicity in this section, we refer to it as problem (P). 

Note: non-positive ;r's and v's have been converted to non-negative /r's and 
v's. 

In our Benders decomposition scheme, the interdiction vector 8 plays the 
role of the vector of “complicating variables.” 

Rewriting (P) as: 

(P): max max c v v + c K 7t 

Se A v,7t 

s.t. A v v + A 71 n = b - BS ( y ) ^2 16) 

v > 0 
n > 0 

the dual of the inner maximization above is: 

min y(b-B8 ) 

s.t. yA' >c v (v) (2.17) 

yA*>c* (7z) 

where y represents the vector of dual variables corresponding to the constraints 
of the primal problem in (2.16). We assume the polytope 
H =\y\yA v >c v ,yA n >c*j always contains a feasible solution, and its optimal 
solution is not unbounded. We will show later that these assumptions are valid. 

Since the solution to (2.17) will lie at one of the extreme points of the 
feasible region H , we can rewrite (2.17) as: 
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min y,(b- BS) (2.18) 

iel 

where / = {1,2,..., A} indexes the set of extreme points of H which is 

{vi,j 2 ,...,v |7| }. 


Since min y.(b-B5) < yXb-BS), Vie I , we can rewrite (P) as: 

iel 


(P): max z 

5eA,zeR 


(2.19) 


s.t. z < y t (b - BS) Vz e / 

If we relax (2.19) by considering a subset of constraints iel, indexing 
these constraints by i = \,2,...,k, at a given iteration k of our Benders 
Decomposition Algorithm (BDA), we have the following “master problem” (MP) to 
solve: 


(MP,): z k = max z 

5eA,zeR 

s.t. z < y^b-BS) 


i = 1 , 2 , 


k 


( 2 . 20 ) 


(MP*) is a relaxation of (P), and therefore its optimal objective function 


value, z k , constitutes an upper bound on z* , the optimal objective value to (P). 


If the solution to (MP t ), denoted ( 8 k ,z k ), is not optimal to (P), then (S k ,z k ) 
must violate some constraints in (2.19) not included in (MP* ), i.e., constraints 
associated with extreme points {y k+K y k+2 ,-~,y , } ■ (Remark: Here, the subindex k 

refers to the iteration counter, and S k is the incumbent interdiction vector.) 

A natural way to check for one of these extreme points is to identify the 
vector y k +\ satisfying y k+x e H that violates z k < y k+l (b-BS k ) the most. That 
identification problem can be stated as the following optimization “subproblem:” 

(SP*(4)): min y k+1 (b-BS k ) 

Jt+i 

s.t. y k+l A v > c v (V*) (2.21) 

y k+ ^>c* (K k ) 


or, in its primal version: 
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( 2 . 22 ) 


(SP,(4» : max c v v + c n 7t 

V,7T 

s.t. A v v+A*7t = b- BS k (y k+ \) 
v > 0 
7l> 0 

It is interesting to note that model (SP*) in (2.22) is the same problem as 
problem (P) in (2.16) when the value of the decision variable 8 is fixed as 8 = 8 k . 
(SP A ) is therefore feasible, because problem (P) is feasible for any 8 = 8 k . 

(Recall that the inner problem in (2.16) is equivalent to a DCOPF model, which is 
always feasible, even if this feasibility entails penalties for unmet load.) In 

addition, (SP A ) is bounded for any 8 = 5 k , because our DCOPF models have 
non-negative cost coefficients only. 

The solution to (SPJ, denoted (v t ,fi k ), along with 8 k from model (MP t ), 
form a combined solution (S k ,v k ,7t k ) , which is feasible to our original (P). Its 
objective function value is cv k + c n 7c k , which represents a lower bound on the 
optimal solution to problem (P). 

Solving the subproblem above yields a new extreme point y k+l , whose 
associated cut z< y k+l (b-B8) will be violated by the incumbent solution ( z k ,8 k ) 
unless ( z k ,S k ) is already optimal to (P). The newly generated cut is added to 
(MPJ , becoming (MP i+1 ), and the process is repeated. 

E. DECOMPOSITION ALGORITHM 

The mathematical derivation described in Section D, leads to the following 
step-by-step algorithm: (The &(■) notation refers to the optimal objective 
function value of the problem in the argument.) 

Benders Decomposition Algorithm (BDA): 

Input: Initial solution matrices A v , A* , and B, vector b, and an 

optimality tolerance e>0. 


20 



Output: s -optimal interdiction plan 8*, and associated power flows, 
generation phase angles, and load shedding (jointly denoted as /). 

1. Set the iteration counter k: =0, the lower bound LB:=-°o, and the 
upper bound UB:=+°o. 

2. Solve SP,(4) for ft,Let LB, :=!>{sp,(4)) . 

3. If LB a <LB, then update the lower bound: LB^LB^, and set 

5 8 k , v := y k+1 ■ 

4. If UB - LB < e, STOP; otherwise continue to step 5. 

5. Add the newly generated cut: z < y t+1 (b - B8 k+1 ) to (MP t ). 

6 . Set k:=k + 1. 

7. Solve (MP a ) for 8 k , and z k =^(MP ;t ). 

8 . Update the upper bound: UB := z k , and return to step 2. 

At each iteration, the solution y k+1 to SP(4) gives us a new candidate 
solution ( 8 k ,y k+l ). Whenever the objective value for this candidate solution is 
greater than the previous lower bound, we replace the lower bound with this 
value. This provides a non-decreasing lower bound. The solution to (MP^) is 

monotonically decreasing with each cut added, and therefore we automatically 
update the upper bound at each iteration. This procedures is repeated until the 
bounds converge, or are sufficiently close, which in the worst case will occur after 
generating all possible extreme points of the subproblem, i.e., when k = K. 

To improve the BDA efficiency, we incorporate some enhancements in 
Chapter VI. Figure 4 depicts the original BDA for ease of comparison with the 
flow chart for the enhanced BDA. 
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Benders Decomposition Algorithm (BDA) 


LB:=-oo 
UB:=°o 
k :=0 


Figure 4. 



Benders Decomposition Algorithm flowchart. 
22 












Using the above template, we implement BDA using the full formulation of 
the interdiction problem, with and without system restoration, and test it on 
several benchmark cases. 

The decomposition for the model without system restoration is presented 
in Chapter IV, starting with a simple example to familiarize the reader with the 
process. Chapter V expands upon this model to take into account system 
restoration, and presents computational results for that model for BDA. 

Before describing details of the decomposition process, we validate the 
inner power-flow model in the next chapter. 
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III. VALIDATING THE DC POWER FLOW MODEL 


This chapter validates the DC Power Flow model (DCPF) introduced in 
Chapter II, whose accuracy must be established to ensure the reasonableness of 
our interdiction model, l-DCOPF. We do this using non-interdicted and 
interdicted instances of the IEEE One Area test case (IEEE Reliability Test Data 
1999-1). 

A. INTRODUCTION 

DCPF is a simplification of the full AC Power Flow model (ACPF) (e.g., 
Wood and Wollenberg 1996). ACPF includes reactive power flows (DCPF 
disregards these), voltage magnitudes at the buses (DCPF assumes nominal 
voltages, i.e., 1.00 p.u.), and power losses on the line (DCPF assumes these to 
be zero), among others. 

Overbye et al. (2004) conduct an analysis to validate the use of DCPF in 
place of ACPF for market analyses; they find it adequate. Wood and Wollenberg 
(1996) indicate that DCPF is adequate for security analysis of many systems. 
Although we believe that similar results should hold for interdiction problems, 
further analysis is desired because these problems have additional 
complications. In particular, they involve removing critical components of grids 
and adjusting loads to minimize disruption costs. 

We validate DCPF through comparative analysis with the ACPF provided 
by software in the PowerWorld Simulator (PWS) (PowerWorld 2003). PWS is 
standard software, widely used in the electric power industry. Our assessment is 
carried out by comparing “optimal power flows” (OPFs) obtained using DCOPF 
with those identified by the ACPF. It is important to note that ACPF does not, per 
se, attempt to minimize load shedding (or its cost) as DCOPF does. Flowever, 
for given loads and generation, it provides accurate active and reactive power 
flows on the lines and bus voltages (magnitude and phase), along with other 
system values. We want to assess whether power flows and voltage angles, 
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obtained through DCOPF, are in fact close to those from ACPF when we use our 
optimized generation and load shedding as input data. 

We present the comparison between both approaches assuming no 
interdiction first, followed by other cases where critical components have been 
interdicted. 


B. DCOPF VS. ACPF IN THE ABSENCE OF INTERDICTION 

To enable comparison between the DCOPF and the PWS ACPF, we 
incorporate the IEEE RTS One Area case into PWS through its one-line diagram 
building interface. First, we reconstruct the underlying power grid, as shown in 
Figure 5, using the following symbols to represent the different system 
components: 


Bus 



Line 


Fraction of line 

Fraction of line 

Generator 

capacity in use 

capacity in use 


(if <80%) 

(if >80%) 




Load 


V-Xl/N/ 

AM 




Transformer 


Flow direction 


Circuit breaker 
(close) 


Circuit breaker 
(open) 


For simplicity (and given that generating units are not interdictable in our 
test cases), we aggregate all generating units at a bus as a single generator at 
that bus. Next, we associate grid data with the power grid. The reconstructed 
one-line diagram includes data for bus nominal voltages, line resistances, 
reactances and thermal limits (line capacities), aggregated generating capacities, 
etc. Then, we incorporate active power output for each generating unit and 
active load met at each bus from DCOPF as system data for PWS’s ACPF. 


26 





Finally, we run ACPF to compute actual power flows, phase angles and adjusted 
generation at the swing bus. The resulting power flows on the lines are 
compared to those provided by DCOPF. 


The one-line diagram in Figure 5 represents a case without interdiction, 
which we call “scenario 0,” because it is equivalent to a case with zero 
interdiction resource. Section C analyzes cases with interdiction. 



Figure 5. IEEE RTS One Area case one-line diagram without interdiction. This 

representation is constructed using the PWS software (PowerWorld 2003). 
It depicts flow directions, fractions of line capacities used, generation 

outputs and loads. 

A side-by-side comparison of the resulting flows is shown in Table 1. We 
analyze absolute deviations for all lines and transformers. The column headed 
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“Line” displays the line name. The “From Bus” and “To Bus” columns indicate 
the two buses connected by the line. The “Transformer” column indicates if we 
are using a line to represent a transformer. The “PWS flow” and “PWS loss” 
columns show active power flows and losses, respectively, obtained with PWS, 
while the “DCOPF flow” column displays active power flows obtained with 
DCOPF. The last column shows the percentage absolute deviation between the 
PWS and DCOPF flows. The results show a maximum deviation of 5.29%. 
Mean deviation is 1.18% with a percentage standard deviation of 0.012 %. 

The example above indicates that DCOPF yields a good approximation for 
the true AC power flow in this RTS case, in the absence of interdiction. In the 
following section, we analyze what happens when critical elements of the grid are 
removed through interdiction. 


Line 

From 

Bus 

To Bus 

Transformer 

PWS 

flow 

(MW) 

PWS 

loss 

(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A1 

101 

102 

No 

16.70 

0.01 

17.00 

1.80% 

A2 

101 

103 

No 

-26.50 

0.42 

-26.73 

0.87% 


101 

105 

No 

54.00 

0.68 

53.73 

0.50% 

A4 

102 

104 

No 

28.60 

0.29 

28.42 

0.63% 


102 

106 

No 

44.00 

1.01 

43.58 

0.95% 


103 

109 

No 

40.10 

0.53 

39.85 

0.62% 


103 

124 

Yes 

-244.30 

1.21 

-246.58 

0.93% 


104 

109 

No 

-45.30 

0.60 

-45.58 

0.62% 


105 

110 

No 

-17.20 

0.07 

-17.27 

0.41% 


106 

110 

No 

-96.80 

0.28 

-92.42 

4.52% 

All 

107 

108 

No 

-123.50 

2.67 

-125.00 

1.21% 


108 

109 

No 

-151.90 

11.66 

-159.94 

5.29% 

HEEV 

108 

110 

No 

-130.50 

8.45 

-136.06 

4.26% 

A14 

109 

111 

Yes 

-153.80 

0.48 

-154.38 

0.38% 

xm . 

109 

112 

Yes 

-185.20 

0.69 

-186.28 

0.58% 


110 

111 

Yes 

-203.10 

0.84 

-204.43 

0.65% 


110 

112 

Yes 

-234.20 

1.11 

-236.33 

0.91% 


111 

113 

No 

-223.90 

3.11 

-226.01 

0.94% 


111 

114 

No 

-132.20 

0.89 

-132.81 

0.46% 


112 

113 

No 

-170.00 

1.78 

-171.01 

0.59% 


112 

123 

No 

-245.20 

7.68 

-251.60 

2.61% 


113 

123 

No 

-183.20 

3.86 

-186.01 

1.53% 


114 

116 

No 

-322.20 

5.36 

-326.81 

1.43% 


115 

116 

No 

45.60 

0.04 

45.90 

0.66% 


115 

121 

No 

-225.00 

3.14 

-227.24 

1.00% 
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Line 

From 

Bus 

To Bus 

Transformer 

PWS 

flow 

(MW) 

PWS 

loss 

(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A26 

115 

124 

No 

248.10 

4.33 

246.50 

0.64% 

A27 

116 

117 

No 

-310.20 

2.96 

-312.52 

0.75% 

A28 

116 

119 

No 

87.40 

0.23 

86.81 

0.68% 

A29 

117 

118 

No 

-173.20 

0.61 

-172.69 

0.29% 

A30 

117 

122 

No 

-137.90 

2.78 

-139.82 

1.39% 

A31-1 

118 

121 

No 

-52.30 

0.08 

-52.85 

1.05% 

A32-1 

119 

120 

No 

-47.20 

0.11 

-47.19 

0.02% 

A33-1 

120 

123 

No 

-111.20 

0.38 

-111.19 

0.01% 

A34 

121 

122 

No 

-158.70 

2.35 

-160.18 

0.93% 


Table 1. ACPF versus DCOPF: IEEE RTS One Area case (no interdiction). For 
each line and transformer, we show the PWS’s ACPF power flow and 
losses. We compare power flow absolute deviations from those of our 
DCOPF model. The maximum percentage deviation is 5.29% for the 
transmission line connecting buses 108 and 109. The average deviation 

is 1.18%. 


C. DCOPF VS. ACPF AFTER INTERDICTION 

In order to further validate DCOPF, we also compare flows after 
interdiction. The main effect of interdiction is, normally, load shedding. We 
conduct the analysis for a number of possible scenarios; recall that “scenario” 
corresponds to the amount of interdiction resource available. The flow- 
comparison procedure is repeated for the RTS One Area case for scenarios 
M = 2, 4, and 6, in l-DCOPF shown in Appendix A. 1. 

Although overall results are presented for all cases, we will illustrate the 
one-line diagram and comparison table for scenario M = 6: We open the 
interdicted lines (found for the optimal interdiction plan which corresponds to this 
scenario) by opening the circuit breakers, at the ends of the lines, which are 
represented by green open squares in Figure 6. Then, we run PWS to realize 
that the flow in a number of lines greatly exceeds the lines’ capacities (red circles 
in Figure 6). While a line can exceed its nominal capacity temporarily, the long¬ 
term emergency rating on a line is typically 20% of its nominal rated capacity, 
and a line can handle this excess flow only for 24 hours (IEEE 1999-1). 
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Figure 6. Effects of interdiction to IEEE RTS One Area case. Some lines become 
overloaded, which is unacceptable but for a few hours and up to certain 
limits. Note: All loads are kept at their original values to illustrate the 
infeasibility of the problem after interdiction. 

Having determined that ACPF shows the situation to be infeasible after 
interdiction, we reduce the loads in ACPF to match those provided by DCOPF 
and run PWS. We then carry out a comparison of the flows across the non- 
interdicted elements of the grid with those provided by DCOPF; see Table 2. 

Note that, again, all the DCOPF flows are relatively close to those 
produced by ACPF. The average absolute deviation is 3.67% with a standard 
deviation of 0.0075%. 
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Line 

From Bus 

To Bus 

Status 

Transformer 

PWS 

flow(MW) 

PWS 

loss(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A1 

101 

102 

Closed 

No 

-3.60 

0.41 

-3.23 

10.28% 

A2 

101 

103 

Closed 

No 

68.30 

0.02 

66.01 

3.35% 

A3 

101 

105 

Closed 

No 

136.20 

1.18 

129.22 

5.12% 

A4 

102 

104 

Closed 

No 

104.80 

1.93 

99.95 

4.63% 

A5 

102 

106 

Closed 

No 

92.60 

9.82 

88.82 

4.08% 

A6 

103 

109 

Closed 

No 

76.70 

0.00 

76.52 

0.23% 

A7 

103 

124 

Closed 

Yes 

-10.40 

0.00 

-10.51 

1.05% 

A8 

104 

109 

Closed 

No 

95.70 

0.00 

99.95 

4.25% 

A9 

105 

110 

Closed 

No 

124.90 

0.00 

129.22 

3.34% 

A10 

106 

110 

Closed 

No 

82.80 

0.00 

88.82 

6.78% 

All 

107 

108 

Open 

No 

0.00 

5.00 

0 

0.00% 

A12-1 

108 

109 

Closed 

No 

-33.70 

0.00 

-35.42 

4.86% 

A13-2 

108 

110 

Closed 

No 

-36.80 

0.00 

-39.58 

7.02% 

A14 

109 

111 

Closed 

Yes 

-27.80 

0.00 

-29.6 

6.08% 

A15 

109 

112 

Closed 

Yes 

-3.90 

0.00 

-4.35 

10.34% 

A16 

110 

111 

Closed 

Yes 

-19.80 

0.00 

-20.89 

5.22% 

A17 

110 

112 

Closed 

Yes 

3.90 

0.00 

4.35 

10.34% 

A18 

111 

113 

Open 

No 

0.00 

0.00 

0 

0.00% 

A19 

111 

114 

Closed 

No 

-56.60 

0.00 

-50.49 

10.80% 

A20 

112 

113 

Open 

No 

0.00 

0.00 

0 

0.00% 

A21 

112 

123 

Open 

No 

0.00 

0.00 

0 

0.00% 

A22 

113 

123 

Closed 

No 

-258.70 

7.81 

-265 

2.38% 

A23 

114 

116 

Closed 

No 

-50.40 

1.89 

-50.49 

0.18% 

A24 

115 

116 

Closed 

No 

204.90 

0.00 

204.49 

0.20% 

A25-1 

115 

121 

Open 

No 

0.00 

0.00 

0 

0.00% 

A26 

115 

124 

Closed 

No 

13.00 

0.00 

10.51 

19.15% 

A27 

116 

117 

Open 

No 

0.00 

7.59 

0 

0.00% 

A28 

116 

119 

Closed 

No 

313.50 

1.29 

309 

1.46% 

A29 

117 

118 

Closed 

No 

122.20 

2.29 

123.55 

1.09% 

A30 

117 

122 

Closed 

No 

-122.20 

0.00 

-123.55 

1.09% 

A31-1 

118 

121 

Closed 

No 

-104.50 

0.54 

-104.73 

0.22% 

A32-1 

119 

120 

Closed 

No 

64.30 

0.01 

64 

0.47% 

A33-1 

120 

123 

Open 

No 

0.00 

0.08 

0 

0.00% 

A34 

121 

122 

Closed 

No 

-174.70 

0.00 

-176.45 

0.99% 


Table 2. PWS’s ACPF versus DCOPF: IEEE RTS One Area case, scenario M= 6. 

All columns are as explained in Table 1 with the exception of the “Status” 
column which shows which lines are open. Note that the largest 
percentage deviations occur across lines whose flow is small. 

Results for scenarios 2 and 4, summarized in table 3, are similar. 
Differences in line power flows are negligible either because of their small 
relative (percentage) values, or because of the low power flowing across lines. 
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Scenario 

Statistics for "Absolute Deviation" 

Max. (MW) 

Max. (%) 

Average (MW) 

Average (%) 

Std. dev. 
(MW) 

Std. dev. 
(%) 

0 

8.04 

5.29 

1.70 

1.18 

1.96 

1.22 

2 

2.37 

13.04 

0.44 

0.92 

0.50 

2.28 

4 

10.37 

4.84 

1.24 

0.83 

2.10 

0.96 

6 

6.98 

19.15 

1.95 

3.68 

2.21 

0.96 


Table 3. Overall comparison of power flows provided by DCOPF and PWS. The 
average deviation across the four scenarios is less than 5%. 

Overall, the differences found are probably acceptable in the context of 
solving interdiction problems. This concurs with the generally accepted notion 
that DCPF renders acceptable results in the context of contingency analysis 
(e.g., Wood and Wollenberg 1996). We assume that the differences will remain 
small in all our analyses. 
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IV. BENDERS DECOMPOSITION FOR THE PROBLEM 
WITHOUT SYSTEM RESTORATION 


This chapter describes Benders decomposition applied to the MIP 
interdiction model LDI-DCOPF without system restoration presented in Chapter 
II, Section B. The chapter begins with the application of the BDA described in 
Chapter II, Section D to a three-bus test case with interdictable elements limited 
to lines. The general case is then addressed. 

A. SMALL TEST CASE PROBLEM DERIVATION 
1. Description 

This section illustrates Benders decomposition applied to a small, three- 
bus test grid. This is intended to familiarize the reader with the mathematical 
derivation for the general case that will be described later in this chapter. Even 
for this small example, the mathematical derivation is somewhat tedious. In 
order to simplify the exposition, we restrict the interdictable components to lines 
only. Figure 7 gives the one-line diagram for the example grid. 


Bus 1 Bus 3 



Figure 7. Three-bus (and three-line) example. Buses 1 and 2 are connected to one 
generator each. Bus 3 is connected to a load of L, MW. All the other 

symbols in the figure represent decision variables: P x u ,P 2 are generator 
outputs; 0 { ,0 2 ,d 2 are phase angles at the buses; P i2 ,P' 2 ,P 23 are power 
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flows on the lines. 


2. Mathematical Formulation 

In this example, generating costs are disregarded, so we seek to minimize 
the load shed at bus 3. The problem development starting at the DCOPF level 
follows. 

Decision variables(units): 

P -’, Generator output at bus z, for z = 1,2 (MW) 

Py , Power flow on line (z, j) , for (z, j) = (1,2), (1,3), (2,3) (MW) 
iS 3 , Load shed at bus 3 (MW) 

<9 ( , Phase angle at bus z, for z = 1, 2, 3 (radians) 

Problem data (units): 

, Demand at bus 3 (MW) 

r .., Xj , Line resistance and reactance, respectively, for (z, j) = (1,2), (1,3), (2,3) (L>) 

x 

B tj , Series susceptance of line (z, j). By = 2 2 , 

r y +Xy 

for (i,j) = (1,2), (1,3), (2,3) (l/G) 

P : °, Maximum generating capacity at bus z , for z= 1,2 (MW) 

, Line (z, y) transmission capacity, for (z, y) = (1,2), (1,3), (2,3) (MW) 

(Recall that all our units are converted to p.u. values as described in 
Chapter II, Section A). 

In this small DCOPF example, we seek to minimize load shed at bus 3, 

that is: 
min S 3 
Subject to: 

Power balance constraints: 

Pf-Pn-Pn =0 
Pf+Pn+P* =0 
Pn+P^ + S, =L 3 
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Admittance constraints: 


P n~ B \oA + B vA =0 

P n~ B \A\ + B \Ai = 0 

P 23~ B 23 B 2 +B 23 B 3 =0 

Power Generation bound constraints: 

if <if 
if < f G 

Line capacity constraints: 

p L < p L 

pL > _pL 

1 \2 — 1 12 

P L < P L 

1 13 - 1 13 

pL > _pL 

J 13 - J 13 

P L < P L 

1 23 ~ 1 23 

pL > _pL 
1 23— 1 23 

Upper bound on load-shedding constraint: 

S3<L, 

Variable sign constraints: 

S, >0 
Pf, P 2 g > 0 

if, if, if unrestricted 
f, (9 2 , 0 S unrestricted 

We can now extend this formulation to account for potential interdiction by 
introducing three variables to represent interdiction of the lines. Note that the 
interdiction variables must turn off or turn on the power flow on the lines (i.e., 
open or close the lines, respectively). When f =1, the line (i,j) is interdicted, 

and when f = 0, the line ( i,j ) is left intact. 
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Letting M n = p!. + , where 0 :j bounds the maximum difference 

between phase angles at buses i and j, the interdiction problem for given 
interdictions variables 5^,5^, and S 23 , can be written as follows: (Note: The 
variables n represent the dual variables of the balance and bound constraints.) 
(I-DCOPF) :max min S 

Se A P a ,P L ,S,S 

s.t.: 


pG pL pL 
r \ M2 M3 

= 0 

(*f) 

^+^12+^23 

= 0 

(<) 

^3 + Pli + M 

=2, 

(^3 S ) 

P\2 ~ B\1@\ + B\2@2 

<M,X 

(< 2 ) 

P\2 ~ B\ 2 6\ + B X1 6 2 


04 + ) 

P\2 — + ^13^3 

<M,/‘ 

04 ) 

^13 _ -®13^1 + ^13^3 

> -M n S,\ 

04 + ) 

P 2 3 - B 2 £ 2 + ZTM, 

< m 2 &, 

04 ) 

^23 — B 2 £ 2 + B 22 6 2 

>-M 1s S l „ 

(<3) 

Pi f 

<P,° 

04) 

P 2 g 

<P? 

04 ) 


P L 

1 12 


P L 

1 12 


P L 

1 13 


P L 

1 13 


P L 

r 23 

s/go-4) (<■' 

4 


S3 

Si, «" a ) 


S 3 >0 

P G ,Pf >0 

P\i i P\i ? ^23 unrestricted 

unrestricted 
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where: 


A = {tf, 


' s,\,& L 


12 9 13 ? 


23 


S 12 + 8 U + 8 


< 1 S\L 

23 — ' 


{ 0 , 1 }} 


because we assume one unit of interdiction resource is available. 

Notice that every admittance equation has been split into two inequalities. 
This is done in accordance with equation (2.9) (see Chapter II, Section B), to 
linearize the multivariable product that results from introducing the interdiction 
variables into the model. When = 0, every pair of inequalities enforces 

Py -Bfl + ByOj = 0. When 8^ =1 (line (/, j ) is interdicted), these constraints do 
not bind because M.. is large. However, in this case Pf will still be forced to be 
zero by constraints lf<Pf(l-8j) and If >-Pf(l-8j). 

The optimal interdiction problem consists of maximizing, by choice of 
interdiction variables 8 ‘:, the amount of power that must be shed in the system. 

As shown in Chapter II, Section B, this max-min problem can be easily 
reformulated as a max-max problem by taking the dual of the inner problem, I- 
DCOPF, listed above. This yields: 

(DI-DCOPF): max max ji 3 L, + M r 8' 2 (7r' n - Jtf ) + M r S' 3 (n' v , - k' ]3 ) 

S 71 

+/Jill-<5}, K.T,',:'"’ -) + pt (1 -4- n’i '"') + L,Of) 
s.t. 

4+4 <0 (if) 

4+4 <0 (if) 

4+4° ad <1 (S 3 ) 

~4 + 4 + 4 + 4 + 7t^ ap+ + Jt£ ap ' =0 (if ) 

~4 + 4 + n\ 3 + 4l + n\f ap+ + ji\. f ap =0 (if ) 

~4 + 4 + 4 3 + 4 + 4 Cap+ + Jt4 ap = 0 (if 3 ) 
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—B ] 2 7 T \2 —B 12 7 T l2 —B u 7T u — B 13 7T u 

<3T 

0 

II 

P\ 2^\2 + P\ 2^12 ~ P 23^23 ~ ^ 23^23 

_jn 

<35 

O 

II 

P\ 3^\3 P\ 3^\3 + B 23 ft 23 + B 23 ft 23 

<35 

O 

II 

S^ + S^+S^ 

<1 

7 t S 

unrestricted 

k l ~, 7T LCap ~, tt g , n L 2 ° ad 

<0 

_L + LCap + 

71 , 71 

>0 

cl cL c*Z, 

M2? M3? M3 

e {0,1} 


The resulting model, DI-DCOPF, contains a number of cross-products of 
the form Sn in the objective function, which we can linearize as described in 
Chapter II, Section B. This linearization yields the following MIP: 

(LDI-DCOPF): max L, + M n (vf 2 - v' n ) + M 13 (vf 3 - v' n ) + M 23 (v 23 - v 23 ) 

Se A,v,/r 

+P g 7t g + P g 7t g + P X1 {7t\ Gap - 7T x Gap+ ) + P^ (n[ Gap -7r £ ap+ ) 

-I- P L ( 7T LCap — 'TT LCap+ \—P L (y, LCa P~ _ v LCa P + \ _ P L (y, LCa P _ v LCa P + \ 

' r 2?> VM3 M3 ) M2 V M2 V 12 ) M3 V V 13 V 13 ) 

— P L (yy LCa P‘ — v LCa P + \ T 'Tr Load 
r 23v23 V 23 / MM 


s.t. 



—S 1 __G 

7T X + K x 

<0 

Off) 

. _G 

71 2 + ^2 

<0 

Off) 

7r S 2+7i L 2° ad 

<1 

(M 

-7T X + 7T 2 + ^f 2 + + ;zf 2 Cop+ + 

= 0 

(%) 

+ 7T* + 7T G + 7T x l + 7T x Gap * + 7l\ Gap 

= 0 

(M 

-7T S 2 + 7T° + 7l L 22 + Tt^ + + 7T L 2 Gap 

= 0 

(P 23 ) 

—B n 7T\2 — B l 2 7r n —-5 13 ^ 13 — B n 7t n 

= 0 

m 

B n 7T u + B n 7T l2 — B 237723 — B 22 , 7723 

= 0 

( 0 2 ) 

P\3^\3 P\3^\3 B 23 ft23 B 23 ft 23 

= 0 

(M 

5,2 + 8,3 + 823 

ve V(S,7U) 

<1 



where: 
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<0 


_L _LCap _G „ L , LCap —Load 

7T , 7T , 71 , V , V ,^T 3 

_/. —LCap + L + LCap + \ 

7T , 7T , V , V v >0 

;r v unrestricted 


eZ cL cl 
C7i 2? C> 2 


€{ 0 , 1 } 


The linearization of constraints ve V(S,7r) is specified in detail in Appendix 


C, but as an example, consider the constraint v L =7r L+ S L , where n L * > 0. The 

y y y y 


linearization uses the following constraints: 


U s =L + c 'L 

V.. < 71- O 

y y y 

l + ^ ^L + 

V.. <7T ij 


f ><-<( 1-S{) 


Trf >0 



The procedure delineated above adheres to the formulation developed by 
Salmeron et al. (2004). This formulation allows us to use any general-purpose 
MIP solution technique to solve LDI-DCOPF, and therefore, l-DCOPF. The 
remainder of this section is devoted to demonstrating the application of Benders 
decomposition to LDI-DCOPF in our three-bus network, following the steps in 
Chapter II, Section D. 


We first construct the subproblem for a given vector of complicating 
variables S k at iteration k\ 


(SP/- ( d k )). max 7T k L k + M n (v 12 v 12 ) + M u (v 13 v 13 ) + Af 23 (v 23 v 23 ) 

V,7T 

+P x G 7T° + P G 7T G + P x L 2 (7T^ ap - 7T^ ap+ ) + P^ (TT^ ~ 7T^ P+ ) 

_1_P L (Tr LCa P' _ TT LCap * \ — P L (m LC “P —T, LCa P \—P L (M LCa P —M LCa P + \ 
" ri 23 V*'23 n li ' r 12 VM2 1 12 ) M3 V*13 M3 > 

— P L (v LCa P —v LCa P*\ 4 - T Tr Load 
/ 23V23 M>3 l" 1 " '‘M'M 
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s.t. 








-S l -rr G 

7ty + ^Tj 





< 

0 

(^f) 

—S | _G 

71 2 + 71^ 





< 

0 

(^ G ) 

—S . _Load 

7U i +7T 3 




< 

1 

(5 3 ) 

-7t 1 + ^2 

+ 4 

+ 4 + 

i -rj-LCap* 

-r/J.J 2 

i ji-LCap 

" l "' t 12 

= 

0 

(^) 

_s , _s 

—^Tj + ^T 3 

+ <3 

+ 4 + 

i sr-LCap* 

+ 7T U 

, —LCap- 

" r 'M3 

= 

0 

(^) 

_s , _s 

-^ 2 +^3 

+ 4 

+ ^23 

i jpLCap* 

’ TJL 23 

i Tj-LCap 

TJL 2 3 

= 

0 


—B n 7T n 


n-B 

13-^13 _ ^13^13 

= 

- 0 

(^l) 

B n 7t n + B l2 7T l2 

~ ^23 

^23 ~ ^23^23 

= 

- 0 

( 6 i) 

5 ,3<3 + fi , 3 <3 

+ 5 23 

, n 23 ^^23^23 

= 

- 0 

(^) 


ve V(d k ,7T) 

Constraints veV(S k ,7r) are linearized as in Appendix C with 8 replaced 

by 4 ■ 


The master problem at iteration k becomes: 


(MP /f ): max z 


<5eA,ZeM 

S.t. 


Z< 8\ 2 Y\2\,k' ^12 (1 ^\2^Tl23, "*"^13 ^nTni,k' ~^\3 (t ^13^^33,k' ^23 ^23^231,i' 

~^23 (I — ^23) 7*233,4' _ •^"12 3\2Yl2\,k' "*"^12 _ ^12 )'/l23,k' ~ ^\3 ^nY\3\,k' "*"^13 (1 _ ^13)^33,A' 
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In the above formulation, the subscript £ for any variable (e.g., ) refers 


to the optimal value of the incumbent variable provided by the solution to 
(SP t ,(40) at iteration k'\k'<k. 


Having developed our master problem and subproblem, we implement the 
BDA in Chapter II, Section D, and test the three-bus case for convergence. We 
use the following values for the test case: 


40 



Z 3 = 5 p.u. 

B tj =5 p.u. 

WA G = 3 p.u. 

P\2^n^23= 4 p.u. 

The algorithm converges successfully to the same optimal interdiction 
plan: (S n = 0,£ 13 = 1,<? 23 =0) as that of the MIP formulation in LDI-DCOPF. The 

following graph illustrates how upper and lower bounds converge after four 
iterations of the algorithm. 


BOUNDS VERSUS ITERATION 



Figure 8. Convergence of Benders decomposition for the three-bus case. The 
optimal interdiction plan, which consists of interdicting line (1,3) in Figure 

7, sheds a load of one p.u.. 

This section has illustrated Benders decomposition using a small example. 
This foundation should help the reader understand the generalization of the 
procedure to any power grid in which buses, generators, substations and lines 
can be interdicted. The extension to incorporate the effect of system restoration 
over time will be addressed in Chapter V. 
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B. BENDERS DECOMPOSITION: GENERAL INTERDICTION PROBLEM 

WITHOUT SYSTEM RESTORATION 

The general interdiction problem should be scalable to let us analyze an 
arbitrary power system. We will not discuss the derivation of the BDA for I- 
DCOPF as done in the previous section, however. Appendix A provides a 
complete derivation of the Benders decomposition procedure for l-DCOPF. The 
derivation is broken down into these steps: 

A.1 shows the interdiction model, l-DCOPF. 

A.2 shows the dual of the interdiction model, DI-DCOPF. 

A.3. shows the dual of the interdiction model after linearization of 8n 
products, LDI-DCOPF. 

A.4 shows the BDA master problem. 

The subproblem at iteration k is LDI-DCOPF with all S’s replaced by 
fixed S k ’s (given by the user at iteration k = 0 and by {MP k ) at iteration k > 0), 

where 8 k e A,\/k . 

Once the subproblem is solved for a given interdiction plan, we retrieve 
the dual values for each subproblem constraint and use them as coefficients for 
the master problem. The master problem at iteration k is displayed in Appendix 
A.4. 

C. RESULTS 

The BDA is implemented in GAMS (GAMS 2003) and solved with CPLEX 
version 8.1 (GAMS-CPLEX 2003) as the underlying solver, and run on a Pentium 
4 laptop computer at 3.0 GHz and with 512 Mbytes RAM. We turn off the 
“presolve” option in CPLEX because it leads us to erroneous duals from the 
subproblem. 

Figure 9 shows the sequence of lower and upper bounds generated by 
BDA for the IEEE RTS One Area case (1999-1), with interdiction resource M = 6 
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and the assumptions made in Salmeron et al. (2003). The optimal interdiction 
plan for this case matches the solution obtained by solving the MIP LDI-DCOPF 
directly. 


BOUNDS VERSUS ITERATION 



Figure 9. IEEE RTS One Area case without system restoration: Convergence. Note 
the significant improvement of the lower and upper bound at early 
iterations. As the number of iterations increases the effectiveness of the 
cuts (given by the upper bound from the master problem) decreases. 
Initial upper bound is 140,000, but the vertical axis on the graph is 
truncated for display purposes. 

The results for this problem show that, by iteration 375, the lower bound is 
already within 7.6% of the optimal solution, while the upper bound is within 12%. 
(For this comparison we use the optimal solution obtained with LDI-DCOPF.) We 
achieve these bounds in approximately 10 minutes; however, it takes 1 hour and 
20 minutes to reach the optimal solution which is 14,048 $/hr. This result shows 
Benders decomposition may obtain sensible bounds in a relatively small number 
of iterations, and an acceptable time. However, reaching optimality takes 
considerable effort. This difficulty will be addressed in Chapter VI. 


43 





Figure 10 shows how the time required to solve the master problem and 
subproblem change as the algorithm proceeds. This figure reveals the 
increasing difficulty of the master problem as more cuts are added with 
subsequent iterations. While the time to solve each succeeding master problem 
of the decomposition process significantly increases, the relative efficiency of the 
new cuts are less significant as shown by Figure 9. 
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Figure 10. IEEE RTS One Area case without system restoration: Time versus 
Iteration . The time to solve each master problem in the algorithm 
significantly increases with subsequent iteration. The maximum time to 
solve any subproblem is less than 0.2 seconds while master-problem 
solution times are as high as 18 seconds. 
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V. BENDERS DECOMPOSITION FOR THE PROBLEM WITH 

SYSTEM RESTORATION 


This chapter extends Benders decomposition to the interdiction problem 
with system restoration. We describe the models involved in the new 
decomposition, and show preliminary results. (Chapter VI will show improved 
results through some refinements in BDA.) 

A. FORMULATION 

This section illustrates the application of Benders decomposition to the 
interdiction problem with system restoration, l-DCOPF-R. The derivation of I- 
DCOPF-R is included in Appendices B.1 and B.2. We will not discuss the 
complete derivation of the Benders decomposition procedure for the system 
restoration case, which is included in Appendix B. Specifically: 

B.1 shows time-period constructs. 

B.2 shows the interdiction model l-DCOPF-R. 

B.3 shows DI-DCOPF-R, which is the interdiction model with the inner 
power flow model replaced by its dual. 

B.4 shows DI-DCOPF-R after linearization of cross-products. We denote 
this model LDI-DCOPF-R. 

Finally, B.5 shows the BDA master problem. 

We start the description after the linearization of the 8k cross-products. 
That is, assuming we have LDI-DCOPF-R as the MIP to which we will apply 
Benders decomposition. 

Again, the subproblem at iteration k is the LDI-DCOPF-R with all 8’s 
replaced by fixed 8 k ’s (given by the user at iteration k = 0 and by the (I MP k ) at 

iteration k> 0), where 8 k e A,\/k, and the dual values for each subproblem 
constraint are used as coefficients for the master problem. 
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The BDA master problem is shown in Appendix B.5. At each iteration k, a 
new cut is added to the master problem. This provides a (monotonically 
decreasing) upper bound z k , and an interdiction plan that will be used in the 

subproblem (SP /c (4» to calculate a lower bound. 

B. RESULTS 

Figure 11 shows the convergence of BDA for the IEEE RTS One Area 
case with system restoration. Note the significant improvement in the bounds 
during the early iterations of the decomposition process when compared to the 
last iterations. 


BOUNDS VERSUS ITERATION 



Iteration 


Figure 11. IEEE RTS One Area case with system restoration: Convergence. The 
upper and lower bounds of this problem converge quickly in early 
iterations. Convergence is reached in 120 iterations and takes 33 
seconds. Also note that the optimal objective value (circled in graph) is 
reached in early iterations, but it takes the algorithm a relatively long time 
to prove it has obtained an optimal solution. 
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Although the case with system restoration is intuitively more difficult to 
solve than the case without system restoration (because of the number of time 
periods, and the related number of additional variables), the former actually 
solves faster in the IEEE RTS One Area case. This could be attributed to the 
fact that, for the problem with restoration, some candidate components, such as 
those with the largest repair times, are more obviously attractive than others to 
become part of the optimal solution. This might allow the BDA to target solutions 
that include those components at early stages, which in turn provides accurate 
bounds sooner than in the case without system restoration (in which the 
attractiveness of all components is more balanced). Our results support this 
conjecture. While it takes over 700 iterations and 1 hour and 20 minutes to solve 
the IEEE RTS One Area case without system restoration problem, the problem 
with system restoration solves in only 120 iterations and 33 seconds. We find 
similar behavior in the RTS Two Area case. 

Figure 12 shows that master problem time by iteration remains almost the 
same for the IEEE RTS One Area case with system restoration. We will show 
that this is not the case when we apply the BDA to the IEEE RTS Two Area case. 


TIME BY ITERATION 



Figure 12. IEEE RTS One Area case with system restoration: Cumulative time. Time 
is shown for the subproblem and for the master problem separately. Note 
that the time to solve the master problem appears to remain stable by 
iteration, despite the increasing number of constraints. 
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To illustrate the potential difficulty in solving the master problem as the 
number of iterations increases, we analyze the results for interdiction scenario 
M = 12 for the IEEE RTS Two Area case. Figure 13 shows the results for a 


2,000-iteration run of this problem. 
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Figure 13. IEEE RTS Two Area case with system restoration: Solution trajectory. In 
2,000 iterations, taking 75 hours, a gap of 16% is achieved. The top left 
graph shows that a near-optimal solution is found at the early stages, but it 
takes many iterations to prove it. The top right graph shows the 
cumulative solution time, by iteration, for the master problem and 
subproblem. Note how the master-problem solution times increase 
substantially as the algorithm proceeds. The bottom left plot shows how 
the optimality gap changes with time. Note the large decrease in the gap 
in the early iterations. The bottom right picture is an expanded view of the 
first 30 minutes of the “GAP VERSUS TIME” plot. The gap is reduced to 
50% in the first 30 minutes. The problem’s lower bound is, in actuality, 
only 2% from the optimal solution value at this time, although the algorithm 
cannot prove this within the first 75 hours. 


Clearly, BDA is impractically slow for this problem. This suggests the 
need for strategies to accelerate BDA’s convergence. The next chapter 
describes and demonstrates techniques we have developed to do this. 
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VI. ALGORITHM REFINEMENTS AND RESULTS 


As seen in Chapter V, long master-problem solution times are a major 
factor in the overall efficiency (or inefficiency) of BDA. This prompts us to focus 
on the master problem and reduce its solution times to accelerate overall 
convergence of BDA. 

A. INTRODUCTION 

This chapter explores the following techniques to reduced master-problem 
solution times: 

1. Solving a relaxed master problem in some iterations rather than the 
standard master problem, 

2. Dropping certain Benders cuts as iterations proceed, and 

3. Not solving the master problem to optimality in all iterations (Sub- 
optimal integer solutions). 

The modified BDA which includes the above strategies 1, 2 and 3 is 
shown in Figure 14. 

The various techniques and associated results are discussed in Sections 
B through D below. A combined strategy is reviewed in Section E. 
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Figure 14. Benders Decomposition Algorithm (with refinements) flowchart. This 
diagram depicts the steps of our enhanced BDA, incorporating the 
proposed strategies to speed up convergence. RMP refers to the linear 
relaxation of the BDA master problem. The flowchart includes the 
refinement techniques that are explained in Sections B through D of this 

Chapter. 

50 














B. RELAXING THE MASTER PROBLEM 

The first technique we try solves a relaxed master problem (RMP A ) in most 
iterations k, and only solves the true (mixed-integer) master problem periodically. 
A relaxed master problem should be much faster to solve, and its use may lead 
to an overall reduction in solution time. We refer to this technique as the 
“relaxed-MP strategy.” 

RMPa is formed by treating any discrete variables as continuous (while the 
rest of the problem remains the same). The optimal value of a relaxed 
maximizing model yields an upper bound on the optimal value of the full model. 
However, the solution to the relaxed master problem is not a feasible interdiction 
plan (unless it happens to be an integer solution). For this reason, the relaxed 
master problem solution cannot be used in the subproblem to obtain a valid lower 
bound; however, the Benders cut generated by the subproblem at that solution is 
valid. In order to improve the lower bound, we must solve the actual master 
problem, and we do this at a specified interval (for example, every ten iterations). 
Unless specified otherwise, in the test cases that follow, our relaxed-MP strategy 
consists of solving the true master problem once every ten iterations and 
otherwise solving RMPa. 

To show the benefits of using the relaxed-MP strategy, we first consider 
BDA applied to the IEEE RTS One Area case without system restoration. The 
original algorithm’s solution trajectory for this problem is presented graphically in 
Chapter IV, Figure 9. Recall that it originally takes about 1 hour and 20 minutes 
to solve this case. The relaxed-MP strategy solves the problem in only 26 
minutes, a 67% improvement over the original algorithm. Figure 15 shows how 
the problem converges. 

The reduction in overall solution time cannot be attributed solely to faster 
solutions of the master problem, however, since the modified algorithm solves in 
only 281 iterations compared to the original 768. It appears that cuts generated 
from non-integer solutions may be more effective than standard cuts generated 
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at integer solutions. We are unsure why this occurs in this case, and in fact, we 
will see later that this behavior cannot be generalized to all our cases. 



Figure 15. Relaxed master problem strategy: Convergence (I). IEEE RTS One Area 
case without system restoration. We solve the full master problem every 
10 th iteration and otherwise solve that problem’s continuous relaxation. 
This problem solves 67% faster using this technique compared to the 
original BDA which solves the (mixed-integer) master problem in each 

iteration. 

Figure 16 shows the solution times by iteration, for the master problem 
and for the subproblem. 
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Figure 16. 
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Relaxed master problem strategy: Time versus iteration (I). IEEE RTS 
One Area case without system restoration. We depict computational 
times when using the relaxed-MP strategy. The top figure shows the 
subproblem solution time by iteration (which averages under 0.15 
seconds). The full master problem is solved every 10 th iteration (these 
are easily distinguished in the bottom plot due to their longer solution 
times). The average time required to complete each master problem 
iteration is reduced to less than one second in contrast to the original BDA 
whose iteration average is 9.5 seconds. 


Figure 16 shows that subproblem solution times are approximately the 
same as in the original BDA (see Chapter IV, Figure 10), whereas the master 
problem solution times have been reduced considerably. The apparently 
exponential increase in master-problem solution times (exhibited by the original 
BDA as iterations progressed) does not appear now: The new master-problem 
solution times appear to increase only linearly by iteration. The average time per 
iteration using the relaxed-MP strategy is 0.15 seconds, much less than the 
average of 9.5 seconds per iteration for the original BDA. 


The results above are for the problem without system restoration. Figure 
17 depicts the solution trajectory for the IEEE RTS One Area case with system 
restoration, using the relaxed-MP strategy. 
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Figure 17. Relaxed master problem strategy: Convergence (II). IEEE RTS One Area 
case with system restoration. The algorithm closes the gap rapidly in 
early iterations but takes many iterations to prove optimality. Detailed 
solving time from the algorithm, not shown, reveals that a 10% gap is 
reached in 12 seconds, but optimality is not reached until 60 seconds. 

Note, however, that this problem takes longer to reach optimality using the 
relaxed-MP strategy than the original BDA. This fact can be attributed to the 
number of relaxed master problems (nine in this run) that are solved before we 
solve the standard master problem in order to update the lower bound. If we 
reduce the interval for solving the standard master problem, for example, every 
3 rd iteration, we would be able to solve this problem to optimality in only 23 
seconds, a reduction of 8 seconds (25%) over the original BDA. 

Figure 18 shows the cumulative solution time for the master problem for 
both the original BDA and the relaxed-MP strategy. Note the reduced solution 
time for the master problem when compared to the original BDA. 
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Figure 18. Relaxed master problem strategy: Time versus iteration (II). IEEE RTS 
One Area case with system restoration. We show cumulative time to 
solve the master problem. The average time per iteration for the master 
problem in the relaxed-MP algorithm is 0.03 seconds (including iterations 
where the master problem is solved exactly), a substantial improvement 
from the original BDA, which is 0.15 seconds. Cumulative time is reduced 

accordingly. 

The advantages of the relaxed-MP algorithm are more significant as 
problems become more dificult. This algorithm helps close the optimality gap 
faster because (1) every relaxed Benders cut is valid and can therefore 
contribute to the reduction of the upper bound, and (2) the relaxed master 
problem is so much faster to solve. Figure 19 shows the effect of this algorithmic 
strategy for the IEEE RTS Two Area case with system restoration. 
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Figure 19. Relaxed master problem strategy: Convergence and Time versus 

Iteration (III). IEEE RTS Two Area case with system restoration. The top 
graph illustrates that the relaxed-MP algorithm closes the optimality gap 
more in earlier iterations. The bottom graph illustrates the master 
problem’s reduced cumulative solution time when we use the relaxed-MP 
algorithm. Note that a logarithmic scale on the time axis is used for clarity. 

In the IEEE RTS Two Area case we can appreciate the impact of the 
relaxed-MP strategy on overall solution time. The original BDA takes over 75 
hours to complete 2,000 iterations, reaching a 16% gap. By solving the relaxed 
master problem for nine of ten iterations, and solving the standard master 
problem every 10 th iteration, we reduce the total solution time to a mere 26 
minutes for 2,000 iterations, and reach a 19% gap. The original algorithm 
reaches a 20% gap at iteration 1,683 in approximately 37 hours; in contrast, the 
relaxed-MP algorithm achieves the same gap in 1,762 iterations but only in 20 
minutes. The average time per iteration for this test case using the original 
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algorithm is two minutes and 15 seconds; using the relaxed-MP algorithm, we 
reduce the average time per iteration to 0.78 seconds. The original BDA, of 
course, takes much more time to achieve such gap. 


C. CUT-DROPPING STRATEGIES 

The relaxed-MP algorithm helps solve certain problems, but larger 
problems need additional techniques like “cut-dropping” if they are to be solved 
efficiently. The goal of cut-dropping is to limit the number of Benders cuts in the 
problem so that the master problem remains efficiently solvable. Of course, 
convergence of BDA may be lost unless care is taken. 

We explore the following cut-dropping techniques: 

- Drop (delete) all “sufficiently slack cuts,” 

- Drop all “non-LP-active cuts,” 

- Dropping the first slack cut, 

- Keep a minimum number of cuts plus all active cuts, and 

- Keep the «-most active cuts. 

The following subsections explain and demonstrate each technique in 

detail. 

1. Dropping All Non-active Cuts 

This strategy attempts to keep only binding cuts (the cuts that remain 
“active” between iterations) in the master problem. In a linear problem, an 
inequality constraint is deemed “active” (or binding) at a given feasible solution if 
the exact equality holds at that solution. A non-active constraint is said to be 
slack. The concept of slackness is important here because though basic 
sensitivity analysis one can show that, at the optimal solution to a linear problem, 
slack constraints can be dropped without changing the optimal solution to the 
problem. Moreover, it is well known that, for linear problems solved through 
Benders decomposition, one may eliminate non-active cuts in the master 
problem at a given iteration without losing convergence to an optimal solution. 
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For a mixed-integer master problem, the above property is also true, but 
checking for active cuts is a more complicated task because a cut might be (in 
fact, will be in most cases) slack at the optimal integer solution, yet removing it 
could cause a (strict) relaxation of the master problem. Checking for active cuts 
in a MIP requires removing the cut and solving the new MIP (in order to assess if 
the optimal objective function value has changed). Doing this for every cut is not 
likely to yield an efficient algorithm, of course. Instead, our approach consists of 
dropping cuts that are “estimated” to be non-active. We accomplish this by 
removing all cuts at any iteration that have a slack greater than a percentage of 
the optimal value of the current MP. The slack for a cut of the form z<c‘ S 
(where zel, c,£eM' ! ) at a given feasible solution (z,S) is given by 
s = c'-S-z> 0. Then, 5 is compared with the master-problem objective function 
value z. Ideally, we would use (z,S) = (optimal solution to the master 

problem) to evaluate the (relative) cut slackness. However, sometimes we must 
rely on relaxed-MP solutions (see Section B) or even on sub-optimal solutions 
(see Section D) to carry out this comparison. 

We use the RTS One Area case for basic tests and perform several runs 
with various “slackness criteria.” In particular, we eliminate any cut whose slack 
s is greater than 0.1%, 1%, 5%, 10% and 20%, respectively of the current z. 
The results for all these cases are similar and unsatisfactory: It appears that this 
technique eliminates some cuts that are necessary for the convergence of the 
problem. For example, if we use s>0.20z (20%), the master problem bound 
does not improve after the 85 th iteration, impeding this variant of BDA from 
converging (see Figure 20). 
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Figure 20. Dropping non-active cuts: Convergence. IEEE RTS Two Area case 
without system restoration. Lower and Upper Bounds versus number of 
iterations. We only use estimated active cuts in the master problem. 

These cuts cannot close the optimality gap. 

2. Dropping Non-LP Active Cuts 

This strategy estimates the active cuts in the master problem by using the 
active cuts in the relaxed master problem. We keep only LP-active cuts in the 
master problem. We consider a cut to be LP-active if it is active in the relaxed 
master problem. This technique showed no significant improvement over the 
previous cut dropping technique because it also appears to remove too many 
cuts from the master problem. For example, in the IEEE RTS One Area case 
with system restoration, only four LP-active cuts remain in iteration 566, yielding 
a weak upper bound. 

3. Dropping the First Slack Cut 

We check all incumbent cuts at any given iteration until we find one whose 
slack exceeds a pre-specified threshold, for example, being 20% above the 
master problem objective function value. If none of the cuts satisfies this 
criterion, none is dropped. In addition, a cut that is removed is no longer 
considered in the master problem. 
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This strategy still removes cuts that are needed for convergence. In the 
IEEE RTS One Area case with system restoration and a slackness threshold of 
20%, only nine cuts remain in the master problem by iteration 300 and the 
optimality gap is 597%. When we increase the criterion to 80%, 52 cuts remain 
at iteration 300, but the optimality gap is still unacceptable at 180%. 

A variant of this strategy consists of selecting the first cut whose dual 
variable (at the optimal branch-and-bound node for the master problem solution) 
is zero for elimination. Again, this alternative does not provide us with the 
necessary set of cuts to close the optimality gap. 

4. Keeping the w-Most Active Cuts 

The results in the previous subsections show that we must sustain an 
elevated number of cuts in order to avoid destroying convergence. Additionally, 
we can see from those strategies that if we eliminate cuts based on a pre-set 
slackness criterion we may be eliminating cuts needed for convergence. 

Here we try a new strategy where we limit the number of cuts to a pre¬ 
specified value, n. We assume the master problem with n cuts is solvable in a 
reasonable amount of time. At every iteration, we replace the cut with largest 
slack by the cut generated at the incumbent iteration. This strategy guarantees 
that the “best” n- 1 cuts are always used, keeping the problem to a manageable 
size that will solve relatively quickly. 

Testing indicates that being too restrictive in the value given to n can 
prevent the algorithm from converging. For example, Figure 21 shows that for 
the IEEE RTS One Area case with system restoration, restricting the number of 
cuts to n = 30 is not efficient, whereas the algorithm converges successfully for n 
= 110 . 
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Figure 21. Keeping the n most active cuts: Convergence. IEEE RTS One Area case 
with system restoration. Gap versus iteration plot. Limiting the number of 
cuts makes the problem simpler to solve. However, as illustrated in this 
figure, too few cuts will prevent the problem from closing the optimality 

gap. 

Specifically for n = 30 cuts, we can only close the gap to 60% after 180 
iterations (in 45 seconds), and cannot improve after then. For n = 60 cuts, we 
can achieve a 34% gap in 62 iterations (19 seconds). For n = 90 cuts, a gap of 
4% is achieved in 99 iterations (30 seconds). Finally, for any n = 110 cuts or 
more, the problem solves to optimality in 117 iterations (35 seconds). 


D. SUB-OPTIMAL INTEGER SOLUTIONS 

In this section we consider a strategy that is similar to the master-problem 
relaxation described in Section B. This extension involves establishing a 
termination criterion for the (full, mixed-integer) master problem before it is 
solved to optimality. Typical termination criteria are based on limiting any of the 
following: the number of integer solutions found during the branch-and-bound 
process, the number of nodes explored in the branch-and-bound tree, the 
computational time spent, or combinations of the above, such as: “Stop after a 
maximum number of seconds, if at least one integer solution has been found.” 
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If the integer solution found is not integer-optimal (due to the early 
termination according to pre-specified criteria) then the upper bound cannot be 
updated. However, it is still a candidate solution for improving the lower bound, 
which will be obtained after solving the subproblem for that candidate solution. 

To ensure that the upper bound eventually improves, at least every m 
iterations (e.g., m = 50) the full master problem needs to be solved to optimality. 
(Recall that the upper bound can also be updated whenever we solve a relaxed 
master problem). Figure 22, shows the bound trajectory graph for the RTS Two 
Area case where (a) the master problem is solved to optimality every 50 th 
iteration, (b) the master problem is solved to a sub-optimal integer solution every 
10 th iteration that is not a 50 th iteration, and where (c) RMP* is solved otherwise. 
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Figure 22. Master problem sub-optimal solution strategy: Convergence. RTS Two 
Area case without system restoration. This strategy uses relaxed-MP nine 
of every ten iterations, sub-optimal integer solutions to the master problem 
at every 10 th iteration except every 50 th iteration when the full master 
problem is solved to optimality. After 2,000 iterations, we achieve a 65% 
gap in 26 minutes. In comparison, the relaxed-MP algorithm alone (where 
master problems are solved to optimality every 10 th iteration) achieves a 
62% gap in 1 hour and 42 minutes in the same number of iterations. 
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This strategy performs better than the original BDA at early 
iterations (gap and time per iteration are reduced). However, it reaches a point 
where the convergence stalls. 

E. COMBINED STRATEGY 

Here, we present the relaxed-MP strategy combined with the cut dropping 
strategy that keeps the n most-active cuts. This combined strategy proves to be 
the best among all possible combinations. We show this combined strategy 
applied to the RTS Two Area case with system restoration. 

Table 4 compares results for 2,000 iterations using the different strategies. 
We take « = 500 most active cuts, and then combine it with the relaxed-MP 
strategy where the full MPa is solved at every 10 th iteration only. 


Technique 

Original BDA 

Relaxed MP 

Best 500 cuts 

Combined 

Final GAP 

16.8% 

19.1% 

48.6% 

8.3% 

CPU time 

75h 12 m 

26m 

1h 45m 

16m 


Table 4. Combined strategies: «-most active cuts with relaxed-MP. RTS Two Area 
case with system restoration. This table shows the Gap and Time 
comparisons for the most effective strategies. The results show over 99% 
improvement in time and a 50% improvement in the gap achieved when 
using the combined strategy in lieu of the original BDA. 

In fact, we can even obtain better solutions for this case by increasing the 
number of cuts allowed in the solution strategy. For example, if we increase the 
number of most active cuts to keep to 700, at the end of 2,000 iterations, a gap of 
only 4% remains, but the time to complete this process increases to 20 minutes. 
We can also increase the number of iterations to decrease the gap. Figure 23 
illustrates results for the combined strategy with n = 700 when applied to the RTS 
Two Area case. 
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BOUNDS VERSUS ITERATION 



Figure 23. Combined strategy, n-most active cuts with relaxed-MP. IEEE RTS Two 
Area case without system restoration. 700-most active cuts are used, and 
the full MP is solved every 10 th iteration. The problem reaches a 5% gap 
within 20 minutes and proves optimality in 45 minutes. 

The solution obtained through this combined strategy represents a 
99.12% improvement over the time required by the original BDA. (Comparison is 
made for a 16% gap, which is the gap achieved by the original BDA.) 
Additionally, this combined technique allowed us to prove optimality through 
complete convergence. 


64 




VII. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER 

RESEARCH 


This chapter summarizes the most important findings of the thesis, and 
proposes areas for future research. 

A. CONCLUSIONS 

This thesis has expanded upon the models and solution methods of 
Salmeron et al. (2003, 2004) for optimal interdiction (by terrorists) of electric 
power grids, using limited offensive resources. 

Our first task was to validate the DC power-flow model that forms the core 
of our interdiction models. We compare power flows computed through the full 
AC power-flow model provided by the PowerWorld Simulator (PowerWorld 2003) 
to those computed by our own DC model, DCOPF. The IEEE One Area 
Reliability Test Case (IEEE 1999-1), with several versions representing different 
interdiction scenarios, make up the set of test problems. The average deviation in 
power flows across all scenarios is less than 5%, and all lines showing deviations 
over 10% carry a negligible fraction of the system’s total power. These results 
indicate that the DC power-flow model is acceptable for interdiction analysis. 

After consolidating the mixed-integer formulations of the interdiction 
models (with and without system restoration) developed by Salmeron et al., we 
have demonstrated that Benders decomposition—this involves iterating between 
a mixed-integer master problem and a linear-programming subproblem—is a 
viable technique for solving these problems. We use interdiction decisions as 
“complicating variables” and develop the decomposition methodology through a 
small example first. We then extend the procedure to a generic power grid, apply 
it to larger test problems, and find that convergence is too slow for practical use. 
For example, an instance of the IEEE Two Area Reliability Test Case (IEEE 
1999-1), which has 48 buses, requires two-thousand iterations to close the gap to 
16%; and that takes 75 hours on a 3.0 GHz personal computer. 
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Not only is the number of iterations large, but each one requires us to 
solve a mixed-integer master problem that becomes more difficult to solve as the 
iterations proceed, i.e., as more Benders cuts are added to it. For example, in 
the problem mentioned above, initial master problems solve in less than one 
second, but by iteration 1,900 they are taking up to up to 600 seconds to solve. 
(Subproblem solution times remains stable and short throughout for this problem 
and all others tested.) 

In order to improve the efficiency of the original decomposition algorithm, 
we propose some refinements to the way the master problem is solved. Among 
the techniques investigated, this combination works best: (a) Solve the full 
master problem only periodically, say, every 10 th iteration, and otherwise solve 
the master-problem’s linear-programming relaxation, and (b) drop certain 
“unimportant” Benders cuts to limit the size of the master problem (we keep at 
most n cuts, those that are estimated in some way to be “the most important.) 
Computational times drop dramatically with this strategy. For example, a 4.7% 
gap is reached in 2,000 iterations for the previously mentioned problem, but this 
is accomplished in a mere 20 minutes. Thus, our improved techniques represent 
important steps toward solving large-scale, real-world interdiction problems. 


B. RECOMMENDATIONS FOR FURTHER RESEARCH 

We identify the following areas of research from which the existing work 
could greatly benefit: 

1. Further reductions of the computational burden imposed by the 
master problem. Further research is required in the areas of 
valid inequalities and cut-dropping techniques, among others. 

2. Enhancements to speed solutions of the Benders subproblems. 
The current implementation does not take advantage of the 
decomposable structure of these problems, which arises when 
the goal of the interdictor is to maximize total unmet demand for 
energy (unmet demand for power, integrated over the time) and 
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when repair times for interdicted components can vary. This 
enhancement is not important for small test problems, but it 
should not be disregarded for larger, real-word problems. 

3. Testing on real-world problems. 

4. Embedding of our bi-level model into a tri-level “system- 
protection model.” The ultimate goal of the study of electric-grid 
interdiction is to help analysts develop plans that minimize the 
potential for system disruption. Formulation of a tri-level 
system-protection model, and provision of solutions through 
exact and heuristic methods are thus required. The system- 
defense model introduced by Israeli (1999) can serve as the 
foundation for such work. 
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APPENDICES 


APPENDIX A: MIP MODEL WITHOUT SYSTEM RESTORATION 

A.1 The Interdiction Model (l-DCOPF) 


This interdiction model, l-DCOPF, attempts to maximize the interdiction 
cost (that is, the sum of generating costs plus load-shedding penalties provided 
by DCOPF) by choosing an optimal set of components to be interdicted. The 
following notation is required in addition to that specified for DCOPF: 

Additional sets: 

A*: Interdictable lines (directly or indirectly) 

G *: Interdictable generators (directly or indirectly) 

Additional model data: 


y = ( lif/sr 

[0 otherwise 

i= = { lif s £G ‘ 

0 otherwise 


f 1 if/e ZT 
[0 otherwise 



1 if ge G** 
0 otherwise 


i(g) = Bus for generator g 


s(i) = 


f Substation for bus i, if any 
[0, otherwise 


[ substation for bus i that generator g is connected to, if any 
Remark: s(i(g)) = < 

[O, otherwise 
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I-DCOPF: 


max min 

Se A P Gen ,P Line ,S,6 


ZZ*<0 + Z Z/«<) 

/ gG G t i c 


s.t. 


Admittance constraints, prior to linearization: 

n d-^*) n n <«-#*•) « 

ie/'|teif“ ieS*|feif 4 lleL’\IleL? ar 

which are re-written below after linearization of the SO products: 


p l u --B l (e m -e m )<MM u -s^+ Y. s ‘"+ Z <*”+ Z K‘") vi «) 


iel \IeL, 


seS'\leL s “ b 


lleL I llel?" 


p, u --B,(e, m -e m )>-M,(A l ‘'sr+ Z <T + Z s ?‘+ Z <T‘)vmV) 

Power balance constraints: 


fe/ i/Eir 


//gL |UeLy 


z v - z p . Lu " + z A 1 "+z s *=z^ v/ <*“) 

geG, |]p(/)=i /|rf(Z)=f c c 

Line capacity constraints: 


-pLine ^ pLine 
r i ~ r i 

v/ezT 


pLine ^ ^Line ^ 

V/e / 

{7t\ Cap ' \ 

nLine ^ fiLines-t s?Bus\ 

P, ^ P, (1 -o t ) 

V/,//e/*,/eZf“ 

«) 

pLine ^ pLine ^ £Sub ^ 

\/l,s\seS\leL s ; b 

«) 

rtLine ^ j^Line /i s?Line\ 

p, ^ p i G-s„ ) 

V/,////elV/£ir 

(*£) 

pLine pLine 

r i — r i 

v/e r* 

(^, 4+ ) 

pLine ^ pLine ^ ^Line^ 

V/e / 

(rf Cap+ 

T^Line \ Line /i \ 

P, Z-P, (1-3 ) 

V/, / / e /,/e L] us 

«f) 

pLine ^ pLine ^Sub ^ 

\/l,s\seS\leL s ; b 

Kf) 

pLine ^ "pLine ^Line ^ 

V/,// //g L,lie L p l ar 

« + ) 
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Power generation constraints: 


pGen ^ pGen 
g ~ g 

Vgg G** 

(<°) 

pGen ^ pGen /i cGen \ 

g ~g ^ " °g ' 

Vge G* 

(<) 

pGen pGen /i cBus \ 

g ~ r g °i{g)> 

Vg|/(g)e/* 

(*f) 

pGen < pGen, x _ S Sub j 
g g v *(‘(g)) y 

Vg\s(i(g))*0,s(i(g))G S* 

(*f) 

Upper bounds on the load shedding: 


4 < 4 

Vz,c 

/ _Load 

(% 

Variable sign constraints: 



pGen > q 

Vi,g|ge G t 


pLine URS 

V/ 


O 

Al 

O 

Vz,c 


6. URS 

Vi 



where 

S?Line cBus cSub s?Gen _ f a 1 11 

\8, ,S t ,S s ,S g e {0,1)| 

A < V ' \ ' m Gen <^Gen,t V ' pj Line. I £ Line.t \ ' y j Hus.! £ Bus,t _|_ \ ' Subj <^Sub,t ^ > 

A—i A—i >’g GS Z—l 1 I A—i i i A—i s * — 

G / geG* leH iel* seS* 

A.2 The Dual Interdiction Model, DI-DCOPF 

Taking the dual of the inner model presented in section A.1 yields the dual 
interdiction model shown below. 
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DI-DCOPF: 


max 


Z M < 


W" «-<*)+ X 8?"(zf-zf) 

+ Z «-<)+ Z s “~ «-<*) 

l,s\seS\leL s “ b 7 ,//|//eL*,//eZf“ r 

+Z<" <Z4)+ E/>?+E?(i-‘fK‘+ Z /"o-C)*; 

' e g£G** geG* g|7(g)e/' 


GB 

g 


+ z . 

4 G ( i 

- W 5 ® + V 


-X 

)+Z^ 

g|i(/(g))eS 


/«!** 



leL* 

+ z 

~pLine 

r i 

/i gBus\ /_ LB' —LB + \ . 

(l-A ) O/ j ) + 

z 

Pf ne { 1 

7,i|ie/\/eZf“ 



mb 

+ Z 

r>Line /i s?Line\ /—LL 

P, (1-4 )(%,- 

X« + )+IX 

—Load 

”ic 


l,ll\lleL ,77eZ, i,c 

s.t. 

Power-generation dual constraints: 

*£+d -- 4 Go X° + 4 G <+4g><+4 ( 7(g))< 5 * 4 

Line-flow dual constraints: 

nf + nf + (1 - Af )(nf + nf ) + + Afnf + ^ [nf + nf *) 

l,i\iel* JeLj 

, V (~LS~ . _ LS + \ . V / _Zi7 . _ LL + \ _Bal . _Bal r\ 

+ L, ] + L \ K m +7r Ul j-^(O + ^(/) =0 

7s||eS\7e4 /,//|//ei*,//£if' ,r 

Power-shedding dual constraints: 

_Bal , _Load / 

^ +^i,c ^ he 

Phase-angle dual constraints: 

- Z + ^) + Z + ^ 4+ ) = 0 

l\i=o(l) l\i=d(l ) 

Signs on dual variables: 


Vg 


V/ 


V i,c 


V/ 
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77 BaI unrestricted 


-4 __ Z/n I ■ 

77 ,77^,77 , 

—LB 

n 

-LS 

,77 

_ 7 I 

<0 

77 , 

/rr LB + 

n 

—LS + 

,77 

_zZ 

>0 

_G _G 0 GB 

77 ,77 ,71 . 

-GS 

, ;r 

<0 




n Load < 0 


Bounds on the dual variables: 

n °,7U ,7U ,77 : -77 =-Max. shedding cost-Max. Generating cost 

7r'" , 77 r ,7 T lb ,77 ls ,77 ll : -k l =-Max. shedding cost-Max. Generating cost 

7l L \n L ,7 T lb ,7T ls ,77 ll : 77 1 =Max. shedding cost + Max. Generating cost 

n A : -k a =-P l Line -b,0 

k a : n A =P l Li ' ,e +B,0 

where we take 6 = 1 radian. 

A.3 Linearization of DI-DCOPF 

Linearizing the 8n products in DI-DCOPF yields: 


LDI-DCOPF: 


max^MJ ZOf -vf )+ Z (v BA -v BA+ )+ Z (v" -vff )+ Z 0/“ -v“ + ) 


,/eZf 




/,//|//eL JleL[ 


+Z + Z (< - ) + I (<* - vf) + I P/(< s -vf) 


i|i(i(g))eS 


+x ^ (<° - ^ )+Z “ [of - < + ) - of - v f + )] 

/gzT tezT 

. X 1 j-)Line / —LB —LB + \ / LB LB + \ , X 1 r\Line / —LS —LS + \ / L5 1 LiS + \ 

+ Z 0 O,, )-0 M -v w ) + Z 0 O,, )-0,, s -v /jS ) 

/,s|s£S’,/££f i> 

+ x -<)-(v^-v,“')l+x(l4) <"+ZZ4 

■ ,* ,, ,p„r L_ — I - 


l,ll\lleL JleL[ 


7r + (1 - A 0 )r° + X 7t + A! ji + a rr < h 

i(z) v £ ' £ Z Z i(z) Z s(i(z)) Z Z 


Vg (if) 
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nf + 7if + (1 - A/" 0 )(7rf° + 7rf °) + Afrf + Afxf + ^ (rff + 7t\f j 

l,i\iel* JeL; 



, V (~LS~ . _ LS + \ . V / _ZZ' . _ LL + \ _Bal . _Bal r\ 

+ Z )+ Z +x,,u j-^od) +7r d(l) = 0 

l,s\seS’,leL s ;,;/«eL*,«eif“ r 

V/ 

(tf) 

—Bal , _Load / _/• 

^ +^,c ^ Jic 


V /,c 

(^c) 

- Z ^/(<+<’)+ Z 5 /(^" 

//=o(/) /i=d(0 

+ xf ) = o 

Vi 

W) 

Linearizing constraints are: 




^ ^ gLine 

V t < 7T ; O l 

V/e L* 


(^ + ) 

o 

VI 

+ 

1 

+ 

^ — 
> 

V/e L* 


(/< + ) 
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o 

VI 

+ 

l 
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P^ 
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(^i) 
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v i,s ^ x, 5 S 
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(C) 

o 

VI 

+ 

l 

+ 

^ *5 

P^ 

Vs,l\se S*,l<=L s s ub 


(^ s + , 2 ) 
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o 

VI 
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1 

+ 

3 3 

> 
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07,T) 

IN 

VI 

N 

V/,//|//eZ*,//eZf a '' 

inS) 

IV 

l 

ill 

V/e L* 

inf) 

it 

IV 

i 

in 

ti 

V/e L* 

(Vi) 

—LB - v* =L 

77 H 

V/,/ ie I*,le L Bus 

(C) 

it 

IV 

i 

ii 

\/l,s\seS*,le L s ; b 

(€) 

IN 

1 

Al 

N 

V/,//|//eZ*,//eZf ar 

(vS) 


Vg g G* 

(n G g ) 


Vg|/(g)e/* 

(O 

*?>-K 

Vg s(i(g))e S* 

fof) 

7C Ba> unrestricted 




—A —Ln L _ LS —LL ^ A 

^ ,7T 0 ,71 ,71 ,71 ,71 >0 


v 4 y , v LB y s y L 


>0 


apVl L,, A is is / A 

n ,n a ,n , n ,7i ,tc <0 
, v 1 ", v“~, v LS ~, v LU , v G , v GB , v GS < 0 
7T G ,7t G \7l GB , 7C GS < 0 
7t U,ad < 0 


A.4 The Master Problem 


The master problem for Benders decomposition derives from the dual of 
LDI-DCOPF. At iteration k, the master problem has the form: 

(M P,): max z 

JgA,zgR 
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s.t. 


z < +«,. -KKr-*!%,) 


leL 

l,i\iefjelf* 

+ X (^; 4 ,a- ■ + ^ + )3 ^ nX), ^ > 3 ,a- •) 

l,s\seS\leLf‘ b 

i \ rLine / =A *XA + , =A-J.A + A - J.A -=A *,LA~ \ 

+ X °n ^ n i Yv,u\, k ' +n i Y^uihf-x, J\ unitk <-&I 

I,ll\lleL" JleLf" 

+X K ine (K K,k ■+K Kr - K ?£*■ - K Kx) 

I el' 

+ ^ 8 Bus ( k l y LB+ + jc l y LB — H L y LB — H L y LB 'i 

l,i\iefjeL?™ 

+ ^ S Suh ( 7t l y LS + 7 t l y LS — 7 t l y ls — 7t l y LS ) 

l,s\seS\leL s “ b 

, S Line (jr L y LL +7r L y LL -7r L y LL —7t L y LL X 

/L U ll y yb l 'U,U)„k' yl, l f(Ul) 3 ,k' '*7 !(Ul\,k' i(l,ll) 3 x' 

/,//[//eI* ,lletf ar 

+ z Z 

g\g^G g\geG 


+ Z.C,(-^*'-^C')+ Z/°C' 

g|<(g)eA g|'(g)eA 

+ z z *&& 

g\s(i(g))eS g\sU(g))eS 

+ H<(Kx~Kx > > + X ^(^^-^ 3 ,*')+ X 

l\l ei* /,/|fe/*,/eif“ /,s|.seS\/eZf s 

+ X ) + X (*£*• - ?£*•) + X (?&,*' - ^1,a0 

/,//|//eI* Jletf" l\l eZ* 

+ X (Y(l,s) 3 X ~ Y(l,sh,k'} + X ^7 (Y(l,ll%,k' ~ ¥(1,11%,k’') 

l,s\seS',leLf‘ b 1 ,U\UeL' 

+x^%i-7i))+(^(75--7a0)l+ X -C')) 

+ X [(^(Cu'-Cu'))l+ X [(^(^(W'-C/XA')) 


/,5|5€5 ,/eI, 


/,//|//eZ ,//€i. 


+ Z(-^*x+ z (-«')+ z (-^0+Z«°'+Zx,A 


(i,c),k' 


geG 


i(gW 


s(‘(g))eS 


VA:' = 1,2,...A: 
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APPENDIX B: MIP MODEL WITH SYSTEM RESTORATION 


B.1 Constructing Time Periods, and Additional Notation 

The notation below and time period construct algorithm are extracted from 
Salmeron et al. (2004). 

Required notation: 

T = set of periods, for te T 

^ = i‘uG‘uB*u5‘, set of all (directly) interdictable elements 
Dur(e) = Duration (hours) of outage for element if attacked 

D t = Duration (hours) of time period t, for te T 

1, if component e remains unrepaired in time period t after being attacked 
0, if component e is repaired before time period t after being attacked 
for te T, ee £. 

Remark: In the above notation J3jf, #® us , , and denote J3 te 

when e=l is a line, or e=i is a bus, or e=g is a generator, or e=s is a substation, 
respectively. 



The following algorithm constructs the set of time periods, T, based on the 
different outage durations for all interdictable elements. In the course of the 
algorithm, D t and i) te are also constructed: 
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Algorithm “Construct Time Periods”: 

Initialization: £ <— T 0, t 0, m 0 0; 
While g ^ 0 : 

t i — t + 1 
T^Tvj W 

[0, otherwise 
m t <r- min|z)Mr(e)|ee 

<— jejee £ a Dur(e) = m t J 

End While 


Additional notation for interdiction model: 

C* = Set of lines / that can be directly or indirectly interdicted in period t. 
Note : 

/ e L* if either: 

A!i ine = 1 ’ or 

#* us = 1 for some z|/e or 
/? t s " b = 1 for some s\le L s “ b , or 
j3^ c = 1 for some //1 // e Zf r 
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K 

Kj 

K 


Jl if/e L * 

[o otherwise 
|l if/e lT t 
[0 otherwise 
flifgeG* 

0 otherwise 


*»= HfssC?; ' 
[0 otherwise 

X‘ = j llf,G/ 

0 otherwise 


[l if 5G 5* 

[0 otherwise 


B.2 The Interdiction Model, l-DCOPF-R 


l-DCOPF-R: 


max G min Z D (0-i Z Vf*" + Z Z A A 

s (P Gen ,p me ,s,,e,) t=x [ g t c 


S.t. 


it -w, m -<w * + x p??sr+ z /c<r+ 


is/ /elf 


seS /eZIr 


Z pLine cLine\ 

P,,„ °u ) 


viyt 




lie Lille I?™ 


s e ,o”+ 


iel lie Lf us 


seSlleL s s ub 


Z pLine cLine \ 

Pt,„ 5 „ ) 


V/,Vf 


UeL \lleLY 


IC-Z C'+ I C'+ w >' 

geGj l\o(l)=i l\d(l)=i c c 


«) 
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-pLine ^ pLine 
r t,l ~ r i 

vt,uL; 

«?) 

pLine pLine ^y £Line ^ 

\/t,le L\ fip ne =1 

«/) 

r>Line ^ r^Line /i <?Bus\ 

p tJ ^ p i (!-<4 ) 

Vf,/,i|ie/*, 1&L BUS ,P B “ S = 1 

(*£) 

pLine ^ pLine ^y £Sub ^ 

V^,/,5 5G )S*,/g L s ,j3 B “ b = 1 

0® 

pLine pLine ^y £Line ^ 

Vt,/,// // e L\ll eL Bar , ft™ = 1 

0® 

pLine ^ pLine 
r t,l ~ r i 


«?) 

pLine ^ ^Line^ 

\/t,le L\ P t L j ne =\ 

«/) 

r>Line \ riLine/i rBus\ 

P t , >-P, (1-4 ) 

Vt,l,i\iel*, I e L Bus ,p B ‘ <s = 1 


r>Line \ riLine /1 ?Sub\ 

p a 2-P, (1-4 ) 

Vf,/,j|je5*,/eZ J ,>Sg* = l 

«0 

pLi„ e > _^ ( ' 1 . gLine^ 

gt,i,u\iie l\ug L; a \fc; e = \ 

(*£) 

pGen ✓ pGen 
r t,g - g 

vt, g £ g r 

Kl ) 

pGen ^ pGen /i eGen \ 
r ug — r g V U g ) 

V^,g|ge G\/3°* n =l 

«) 

pGen ^ pGen /-i <?Bus \ 

r t, g - r g y l ~°i( g )) 

Vt,g\i(g)el',/3% g) =l 

«) 

pGen pGen /i <?Sub \ 

r t,g ~ r g ^~°s(i(g))) 

\/t,g\ S (i(g))eS\P™ i(g)) =l 

«) 

S tu * d Uc 

Vz,c 

/_Load \ 

Wt U ) 

pGen > Q 

Vt,g 


P t L ; ne URS 

Vt,/ 


S,J, * 0 

\/t,i,c 


4, OYtt 

Vt,i 



Bounds on : 

7l (,i ' , n G , n GB , n GS : -K ( ' =-Max. shedding cost-Max. Generating cost 
7T l \7T l ,7T lb ,7T ls ,7T ll : -n L =-Max. shedding cost-Max. Generating cost 
n L \n L ,7T lb ,tt ls ,7T ll : n l = Max. shedding cost + Max. Generating cost 
k a : -x 4 =-P l l,ne -B l d 
7T a+ : 7t A = P, Lme + B,e 
where we take 6 = 1 radian. 
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B.3 The Dual Interdiction Model, DI-DCOPF-R 


Taking the dual of the above primal problem yields: 


DI-DCOPF-R: 


i-i-Mw (C-o+ z ft#* «,-<) 


+ z /O” «-<)+ Z MT 

seS’l/eif 6 UeL^etf/" 

+Z(Z4)<'+ZZ^*5+Z Z r°v-sTK, 

i,t c t gid” geG* t\P%r=l 

+ Z. Z Z . Z ?(!-%„)*£ 

g|i(g)e/ AP^“s) =X g|s(«(g))e5 t\p?“\ Kg)) 

+ ZZ^(*5-<i> + Z Z /fa-<f“>«;-4) 

/ /el** /ei* /|/9,j"' e =l 

+ z z p, u (\-s'-)«;:,-<)+ z z jpo-o«,-*£> 

ie/*|/e£f“ f|/?*“=l reS*|/e£f'’ (|/S*f =1 

+ Z Z f7c-a(C-»S)+Zt 

//e£*|/e4'" d/BT=l ',<V 


*&>+<•-« k +s b ( o* ; 


V<,S (O 


*f, +<+(i- A? x*5+ *$ )++ Z (*S + < ) 

!e/*|/e£f“,^®“=l 

+ Z (*£+<£)+ Z (<//+<5/)-<( ; /)+<i) =0 </) 


seS*te£f 6 4 s f =1 


HeLl/letf" 1 


K°+”u?^D(t)-f ic 


V t,i,c (S t i ) 
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- X B I [Ki + 71 a ) + X 5 / (< + </') = 0 

/|o(0=i /|<Z(Z)=i 


VM' (0 U ) 


;r M unrestricted 

_.4 + _Lq _i + _lb + —LS + -LL + ^ A 

K ,K 0 ,K , K ,K ,K > 0 

7Z A ,jf\K L , 7l LB ,7T LS ,7T LL <0 
K ,K ,K , K < 0 
7l Load < 0 

B.4 Linearization of DI-DCOPF-R 


Linearizing the products 8k yields: 
LDI-DCOPF-R: 

n L pLine / A' A + \ . \ 1 pBus / BA' BA + \ 

8/ P,J (v,j - v t l )+ 2^ Pti Gu ~ v o ) 


max 

S. 7 T.V 


X^X 


is/ teZf 


+ X O*- 0 + x /? 7 (C -0 

UeL%L l r! m 


seS \lelz u 


+X(X<)-<"'+XX^ 5 + X X 


z,f c 


‘ ?«<?** 

geG’ i 

1AT=i 

+ X. 

X -o 

X 

X 

g|«(g)eZ 


g|4i(g))eS flA%(*» = 

+XX 

Pfirf-K, 

5 )+X X 

p , L («; 

_I + 
^ t,l 

t i<tC 


/ei* «i^r = 

1 


+ z 

X ^(' 

_ LB + s 

Jt.t,l,i 71 t,!,i ) 

/ LB" 

-0 tj j - 

LB + \) 

V t,I,i )J 






+ X 

X?( 

(-LS- _ZS + \ 

-G tJ ,s ~ 

LS+ \ 
V r,l,s) t 

ssS'l/elf 6 Z|yflf" 4 =l 




+ X 

X 71 

( k ll — k ll 

n t,m 


-v LL+ ' 

IteLllelir t\/5$ e =\ 





r GS -v GS ) 


Load 


4 


With balance constraints: 

^+a- 4 >S + 4 G < +^Xk g X s s * mk 


vt, g (p,d 


86 



K +</ +( 1 -^?)«? +;r ; 

1)+VP,T<, 

, -1 L pLine _L 4 

MA,/ n tji 

+ Z (<u 

i _/-s + \ 







/ 


+ 

Z (<£+<£) 4 

■ Z ( 

—LL . _ LL * \ 

_Zta/ , —Bal A 

-^, 0 (/)+^ ( /) =o 

V t,l 

K) 


/elf 6 ,/?* 6 =1 

//eL*//eZ,f“ r ,/flgf= 1 





_Bal 

x ti ' 

+ ^Z d <D(t)-f ic 




V t,i, 

c{ S tAc ) 

i 

M 

B, (Aj + ^)+Z 5 /( 

K +At)=° 



V t,i 

(*u) 

;|o(/> 

/|</(/) 







Linearizing constraints: 


y4 + ^ ■= A ?Line 

v u 

\/t,i\i&i:,tGT,p ^ e =1 

«/),) 

<;-< 

vr,/|/er,rer,^" e = i 

«/) 2 ) 


\/t,l\l&L\tGT,P^ e = \ 

(</) 3 ) 

V r,/,/ ^,/<7 

Vt,i,l\iel\leLA,teT,P*“ s =\ 


&4 + ^ A 

\/,« ^0 

Vt,i,l\iel\leL? s ,teT,03 a =l 


,5/4 + _yi + s =/4 /i gBus \ 

Vt,i,l\ieI*,leLf a ,teT,fi B J a =l 

(^i )3 ) 

v tj,s ^ x,A 

Vt, S ,l\ S eS*,leL s s ub ,tzT,j3° u s b =l 


SA + *~A + s' a 

Vt,s,l\seS*,leL s ; b ,teT,fi“ b =\ 

(£.o 2 ) 

&4 + _,4 + N =^4 /i gSub\ 

V,j,s -X tJ >~7T tl (\-S s ) 

Vt,s,l\seS\leL s ; b ,t£T,j8°“ b =\ 

(^l) 3 ) 

Z«4 + ^ ^ cLine 

^ x t Ai 

Vt,l,ll\lle L\lle L p ‘ ,r ,t e T,fi£=\ 


LA* „A* s a 

Vt,l,ll\lle L*,lle L Par ,te T,ft ™=1 

(7£i«) 2 ) 

„ LA + —A + \ =^4/1 ?Line\ 

Vf,/,//|//e L*,lle L Par ,tz T,0% e = 1 

(^i« )3 ) 

A / — L ?Line 

v ,,i < K tAi 

Vt,l<=L*,t<=T,j3 b ; ne =1 

(^oi) 

L* „L+ s A 

V r,/ 

Vt,l<=L*,t<=T,j3 b ; ne =l 



Vt,leL\te T,p t L ‘ ne =\ 
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„ LB + / =L gBus 
\i,i 

Vf,/,i |ie I*, Is Lf°,tG T,P?r =l 

(TTO 

LB* —.LB* ^ A 

V t ,l,i <0 

\/t,l,i\iGl*,lGLf s ,tGT,AT =X 

(ATT 

LB + _ LB + \ — L /1 \ 

\/t,l,i\iG f ,/e Zf“ ,tzT,/3*“ s =\ 

(ATT 

I5' + ^ =L cSub 

\,, s < x t ,d s 

Vt,s,l\seS\leL s s ub ,teT,tf?=l 

(&>,) 

v ,.i,s -n tu < o 

Vt,s,l\se S\le =1 

(ATT 

, LS + _Z,.S + \ =1 /i gSub\ 

v <j, s ~X t j, s >-X tJ {\-d s ) 

Vt,s,l\seS*,leL s s ub ,teT,j3™=\ 

(TT h ) 

„ ,ZZ + / =Z, cLine 

v t,i.ll ^ x,jd„ 

\/t,l,ll Ug r,lle Lf ,fe T,AT= 1 

(TinO 

LL + _ZX + / /~\ 

V r,/,// ,1,11 — 0 

Vr, /, // 1 // e L* ,11 e Lf ,t e T, (5'T = 1 

(T,jm )2 ) 

LL + —LL + \ — L /-i cLine \ 

Vr,/,//1 //e L* ,11 e Lf,re T, AT =1 


AT \ ^ cLine 

v u ^-x,A 

Vt,l\leL\tsT,AT =l 

«/>,) 

<;-<> o 

Vt,l\lEL\teT,AT ={ 

(rTJ 


Vt,l\leL\te T,Af e=] 

(rTJ 

TT * -KC 

\/t,i,l\isl\lGLT,tGT,fiT =l 

(TTT 

BA __y4 \ A 

Vt,i,l\i<= I\le Lf“Ve r,#f=i 

(rT ()2 ) 

v ,,u -n tJ <7r lJ {\-d i ) 

\/t,i,l\ie I\Ig L B T ,teT ,#f =x 

(^“•) 3 ) 

^<AT 

Vt, S ,l\seS*,leL s ; b ,t£T,A?= 1 

(&),) 


Vr,j,/|je5*,/eZf i ,rer,^=l 

(& )2 ) 

, ,&4 - _/4 / gSub \ 

v ,j,s ~ K ,,1 ^X u (\-d s ) 

VM,/|seSVeZfVe 7,^=1 

(&) 3 ) 

v tJ.ll ^ -XtAl 

Vf, /, //1 // e L*, I! e Zf’, t e 7, AT = 1 

(TTT 

LA -—A \ /~y 

v r,/,//-^,/ 

Vf, 1,U\UgL\IIg L P T,tGT,AT = 1 

(Y(tJJIh ) 

L^ _ ^ ^ /i oLine\ 

Vf,/,//1 Ug L\Ug zfve t,AT = 1 

(ATT 

vj > 

Vf,/|/er,teT,#f e =l 

(■TuT 

</ -</ 

v?,/|/e r,tGT4:r=i 

(TujT 


Vf, /1 / e Z*, t <e T,/f ne =1 

(TuT 
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LB \ —L qBus 

v ,,u ^-x,A 

VM-,/|ie/*,/eZf“,te 1,^=1 

(T&a) 

v <J,i ~ K tU 

V?,/,/1ZG /*,/e Lf“,te T4 r ; us = l 

(^lz )2 ) 

„ LB~ ~LB~ / =L /i \ 

V?,/,/1ZG /*,/e Lf“,te T,^ r ; us = l 

(AAA 

15“ \ — L gSub 

v u,s ^ -x,A 

Vf,j,/|je5*,/eZf 6 ,rer,^=l 

(AZA 

LS srrLS \ A 

V z,Z,. 

Vf,j,/|je5*,/eZf*,rer,^=l 

(AZA 

, AS *i-LS / =L /i cSub \ 

Vr,j,/|je5*,/eZf 6 ,fer,^=l 

(AZA 

v^z * -KAn™ 

Vf,/,//1 //e L*,ll e L p , ar ,teT, AT = 1 

(AumT 

AZ -<5 ^ 0 

Vf,/,//1 //e L*,lle L Par ,te T,j8 t 'T = 1 

(ATT 


V?,/,// //e L\lle L P T,teT, AT = 1 

(AZA 


Vf,g|ge G*,teT,A e g n = l 

(A, S A 


Vt,g\geG\teT,A; n=l 

(rZA 


Vt,g\geG*,teT,AA=l 

(rf, gh ) 

^-CC> 

Vt,g\i(g)e I*,te T,AZ)= l 

(rZA 

v£-*S*o 

Vt,g\i(g)e I*,te T,AT) =l 

(rZA 

Gi? _Gfi =G /1 rBus \ 

Vt,g\i(g)eI*,teT,AT)= l 

(rZA 

G5 v =G gSub 

v ', s * ~aAsu( S )) 

yt,g\ s (i(g))G s\teT, AT g)) =i 

(rZA 

< g -<:^ 

Vt,g\s(i(g))eS*,teT,A: b Ws)) =l 

(rZA 


Vt,g\ S (Kg))eS*,teT,/}% l(g)) =l 

(rZA 

„A + s ^A 

71 tj ~^t,l 

X/t,l I g Z*,fe T 

A?A 

<z ^ ^Z 

\/t,l\l<=L*,tGT,AT= l 

inf A 

_IB + ^ =Z 
^,Z,Z ~ ^1,1 

Vf, /, 1 1 1 e /V e ZfV e T, Af = 1 

ATA 

—LS + / ;=Z, 

^z,Z, s ^ ^,.Z 

Vf,/„s|seS*,/eZ,f*,fe T,Af= l 

AffA 

_II + ^ =1 
^tJ.U — 71 t,l 

Vf,/,// | //e L*,lle L p , ar ,t&T,AT = 1 

aZu) 
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—.A ^ =A 

71 \,i — 77 i,i 

Vf,/|/e L*,te T 

«) 

in 

i 

Al 

N 

Vt,l\leL\teT,ff; ,e = \ 

(75) 

——LB \ ;=■/, 

71 t,l,i 

\/t,l,i\iGf,lGL Bus ,tGT,/3 B ; s =l 

(C) 

>A Tr 

IN 

1 

Al 

to 

N 

\/t,l,s\sGS\lGL s ; b ,tGT,j3% b =l 

(O 

^LU ^ —L 

71 t,ui - 

1 Ug L\Ug L p ; r ,tG T,p l ™ = 1 


<,>-<* 

Vr,g|ge G*,tG T,ft Ge g " =1 

(O 

a < g 

Vt,g\i{g)Gl\tGT,(3X)= l 

(O 


Vt,g\ s{i(g)) GS*,tGT, = 1 

(O 

n Bal unrestricted 



n A \n^ ,Jt LCap \ 

n LB \n LS \n LL ' >0 



V A \ V ' C “ P \ v LB ,v LS+ ,v LL * > 0 
n r , 7T l ° , ;r iC<,/r , n LB ~, ;r“~, ;r izr < 0 

V A- 5 v IC fl p- 5 v ia- ? v «- ? V LL- ? V G 5 y G5 ? V G5 < q 

TC G ,7T G \7T GB , K GS < 0 

< o 

B.5 The Master Problem 

The master problem of the Benders decomposition results from the dual of 
LDI-DCOPF-R. At iteration k, the master problem has the form: 

(MP^): max z 

JgA,zeIR 
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s.t. 


^ Z Z s i ,ne Kiru,i) v ^ + Kirfw- <r£,i )v k’ -<,ri,) v k) 


t\Pff e =\ leL 


+ X X s r s (Kroii)u 


t\^=U,i\ieI*,leLf us 


+ X X (^/ ^Z), + <1 ^ 1) 3 - </ ^ 1 ), .*' - </ ^ 1) 3 X ) 


^f“*= 1 /,s|.seS* > /elf“ A 


i X 1 X 1 QLine s =A -,LA + . =A ~±A + =A -£A =A ..LA \ 

+ (.^t,iY{t,iji)uk' + ^t,iT(tj,ii) 3 ,k' ^t,iY(t,i,ii\,k' ^t,iY(t,i,ii\,k') 

t\pL™ e =\ IJl\lleL JleL Bar 

r|#r=i 

+ ^ ^ ( y^ L y LB + ft 1 y LB — Ji L y LB — fl L y LB ) 

L-i A-i i v u t,l/(t,l,i\,k' I,l'0,1,>h,k‘ t,ll (t,l,i\,k' yi t,l/ (t,I,i) 3 ,k’> 

t\p b “ s =ll,i\ief ,leL? m 

+ ^ X* § Sub (Jr 1 y LS + jr L y 15 — W L y^ 5 — Jr L y LS ) 

t\p?f=ll,s\s<ES‘MLf b 

+ ^ V 2Line / =L yLL . =L yLL _ =L yLL _ =L yLL \ 

t\tfjF=uji\iteL’ mA" 

+ I Z X’«X M -XX>X + Z Z KX 
+ X X KZ(-<jZ\x-<jZx^ + X X KjZ 

t\P?SzX gl ! '<g) e/ * *\pf“t s x g\‘(s)X 

+ X X ^KC^-C&P + X X <C 3 .* 


G 

r (‘,gh,k’ 


~G n ,GB 

l,g)],k' 


t\fi“U S )XgW<~g^ S 


WT(i( gl X g\s(i(g))£S 


+ Z X^(^)3,r-</)3,r)+ X X K(r^ ih x~K 4 u h x) 


+ X X </(^lv )3 ,r-^l)3,i-')+ X X Kirltuihx-r^uih^ 


t\0™ b =\l,s\seS ,leLf‘ 


t\p% e =\ l,ll\lleL*,lletf ar 


+ X X </0£o3,*'-^) 3 ,*')+ X X 

'IAv”=l W ei * /|^“=l 

+ X X ^t.AY(t,l,s)},k' ~ ^0‘,/,s)3,i-’) + X X ^t,l(7{t,lMh,k' ~ Y(t,l,ll) 3 ,k') 

t\pff=\ I,s\seS‘,leL s s ub APlfl e =1 l ,ll\lle£ Jltl!," 

+ XX^( / 7(i ) ^-<o^) + X llKi^iii ) ,k'- T i(u,\k’) + X X KMtin,k'-vZuk’) 

teT le£ t\pX e =UzL «|y8*"=l /,/|ie/*,tezf“ 

+ X X KM«l)x-vZ,s),k') + X X Ki(€l,nx-€j M k') 

t\p?f=\ l,s\seS\leLf b t\p b f=\ hW^ MI?" 

+ XZK4u') + X X (-<C)P + X X 

4C=lgeG* /|A^)=li(g)e/* 

+ZZ D .K p Lr+ZZ D X s w,w 

teT g teT i,c 


t\p b "1 i(g)) =ls(U,g))eS' 


\/k' = \,2,...,k 
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APPENDIX C: LINEARIZATION OF CROSS-PRODUCTS 


This is an explicit list of constraints for the three-bus example in Chapter 
IV. y,r] represent the dual variables for the linearizing constraints. 


’l2 ^12 ^12 — 0 

(K?i) 

O 

VI 

k 

1 

0&) 

} L — ir L — ir L > —7T L 

12 yi \2 yi \2 u \2 — yi \2 

(^3) 

1 

IA 

0 

04 + i) 

0 

VI 

N 

1 

2 

o&) 


(Ym) 

23 _ ^23 ^23 — ^ 

(Yin) 

£ 

1 

N> hj. 

IA 

0 

(Yin) 

23 ^23 ^23 ^23 — ^23 

iYi D 

’ll +^U,2 >0 

0£.) 

£ “<2 

(r,2 2 ) 

? 12 _ ^12 ^12 ^12 — ^12 

04) 

'l3 + ^13 ^13 — 0 

(^31) 

1 

u> t*J 

IV 

0 

te) 

? 13 ^13 ^13 ^13 — ^13 

(^33) 

23 + ^23 ^23 — 0 

(^31) 

'23-<3 

(^32) 

23 _ ^23 ^23 ^23 — ^23 

(^33) 

LCap + =LCap + ?L < n 

12 '*12 ^12 — v 

oco 

ICap* LCap + / rv 

12 '*12 “ v 

(^ + ) 

LCap + LCap + —LCap + rL > =LCap + 

12 '*12 '*12 ^12 “ '*12 

( Y\iY + ) 
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LCap _ —LCap cL < q 

v u ji 13 u l3 - U 


LCap _ LCap < rv 

V u /£.j 3 i U 


LCap + 
*13 


-n x 


,LCap + 


13 


-X 


LCap 4 


g L > -W LCa P* 

u n - JL n 


(ri 3 


Cap 

131 


( yLCap^ 
^ yLCap^ 


vff* -*ff*<%SO 


V,T- -2-““'’ <0 


LCap 

*23 


LCap + 


-^ _ —LCap _ =LCap ?L > _=LCap 

V 23 Jl 23 Ji 23 U 23 _ ^ 23 


(yLCap + 

(Tv 


Cap 

232 


( yLCap^ 


LCap- , =LCap- CL > q 

V 12 ‘ ^12 ^12 — U 


Off*') 


LCap 

*12 

LCap~ 

*12 


-7Z\ 


LCap 


12 


>0 


_ICap , =LCap cL ^ ^=LCap 

'*12 ^ JL \2 U \2 — ' M2 


(Yl22 


Cap 


(tf 2 


Cap " 
123 


LCap- =LCap- CL > q 

*13 ^ '*13 °3 3 - U 

LCap- _ —LCap- > rv 

*13 '*13 - u 

LCap- _ —LCap- , =LCap- cL ^ = LCap~ 

*13 '*13 " r '*13 ^13 - '*13 


Off* ) 
Off*') 
Off*") 


vff*' +*ff*'4 

>0 


>0 


vff 5 ’' + *ff ,r <5£ £ *ff® 


Of,?*) 


231 

Cap~ 


(T232 


( v LCap 
V/233 


*£ ^ K. 2 

(Jin) 

*S < ng 

(Jin) 

<3 ^ ^23 

(ni) 

<2 ^ -4~ 

(Jin) 

<3 ^ 

(JI 11 ) 

13% 

IV 

A 

(V 23 ) 

—LCap + < =LCap + 

Ji \2 — Ji \2 

(Vn aP+ 

—LCap + / —LCap + 

7U U <7T n 

(Vn^ 

—LCap + < =LCap + 

23 — /<r 23 

(«23 aP+ 
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—LCap _=LCap 
yt 12 - '*12 

—LCap- \ _=LCap 
13 - '*13 

—LCap- > _ =LCap 

JL n — Ji n 


ftf") 

(C') 

(’In’”') 
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